Rasch analysis of outcomes

MADRS-S

Authors
Affiliations
Nils Hentati Isacsson

Centre for Psychiatry Research, Department of Clinical Neuroscience, Karolinska Institutet, & Stockholm Health Care Services, Region Stockholm, Sweden

Magnus Johansson
Published

2024-12-18

1 Setup

Code
library(easyRasch) # devtools::install_github("pgmj/easyRasch")
#library(RISEkbmRasch)
library(grateful)
library(ggrepel)
library(car)
library(kableExtra)
library(tidyverse)
library(eRm)
library(iarm)
library(mirt)
library(psych)
library(psychotree)
library(matrixStats)
library(reshape)
library(knitr)
library(patchwork)
library(formattable) 
library(glue)


### optional libraries
#library(TAM)
#library(skimr)
#library(janitor)
# 935600
# nohup quarto render MADRS.qmd --to html &> render_madrs.out &

### some commands exist in multiple packages, here we define preferred ones that are frequently used
select <- dplyr::select
count <- dplyr::count
recode <- car::recode
rename <- dplyr::rename

path_prim <- '/srv/projects/nils/study4rasch'
path2 <- '/Volumes/projects/k8_CPF_Kaldo/Data/Lärande Maskiner/Nils/Study4data'
running_machine = Sys.info()[['sysname']]
path = ifelse(running_machine!='Linux',path2,path_prim)
set_mywd_path = ifelse(running_machine!='Linux','./study4rasch/rasch_measures','./rasch_measures')
n_cores =ifelse(running_machine!='Linux',1,8)

tryCatch({setwd(set_mywd_path)}, error = function(set_mywd_path){print('already in right directory')})
[1] "already in right directory"
Code
source('settings.R') # containing item labels 
source('mod_easyRasch_func.R') # modified functions + added 

2 Importing data

Code
### import data
df_raw <- read.csv(file=file.path(path,'data','raw_data','outcome.csv')) #

### Clean data 
source("preprocess.R")
df <- preprocess_data(df_raw)

# items to use 
items_to_use <- 'MADRS.S' # options are #MADRS.S, LSAS.SR, PDSS.SR
itemlabels <- itemlabels %>% 
  as_tibble() %>%
  filter(str_detect(itemnr, items_to_use))


# DIF+aux variables to use 
# DIF not used 'marital_status','children', 'working'
dif_variables <- c('Patient','Treatment','sex','age','TreatmentAccessStart','education')

### Make a backup of the dataframe, in case you need to revert changes at some point
df.all <- df

print(itemlabels)
# A tibble: 9 × 2
  itemnr    item                      
  <chr>     <chr>                     
1 MADRS.S_1 Sadness                   
2 MADRS.S_2 Anxiety                   
3 MADRS.S_3 Reduced sleep             
4 MADRS.S_4 Reduced appetite          
5 MADRS.S_5 Concentration difficulties
6 MADRS.S_6 Lassitude                 
7 MADRS.S_7 Inability to feel         
8 MADRS.S_8 Pessimism                 
9 MADRS.S_9 Suicidal thoughts         

3 Checking missing

Code
itemStart <- items_to_use
out = RImissing_mod(df,itemStart,facet="Time")
print(out)

3.1 Remove missing

Code
##### Before filtering out participants, you should check the missing data structure using RImissing() and RImissingP()

# If you want to include participants with missing data, input the minimum number of items responses that a participant should have to be included in the analysis:
min.responses <- 4 # In our case if they have missing on one item all items for that time are missing. 

# Select the variables we will work with, and filter out respondents with missing data
df <- df %>% 
  select(all_of(c(dif_variables,"Time",itemlabels$itemnr))) %>%  # v
  filter(rowSums(is.na(select(.,all_of(itemlabels$itemnr)))) < min.responses) 

4 Create DIF variables

Code
#---- Create DIF variables----
  
# DIF variables into vectors, recoded as factors since DIF functions need this
# these could also be stored in its own dataframe (not a tibble) instead of as vectors
# Named vector for the new types

type_transform <- c(Treatment = "factor", sex = "factor", age = "numeric",
                    TreatmentAccessStart ="integer",education="factor",Time='factor')

# Transform columns based on the named vector
dif <- df %>%
  mutate(across(names(type_transform), ~ switch(type_transform[which(names(type_transform) == cur_column())], 
                                               "integer" = as.integer(.), 
                                               "character" = as.character(.),
                                               "factor" = as.factor(.), 
                                               "numeric" = as.numeric(.)))) %>%
                                               as.data.frame(.) %>%
                                               select(!all_of(itemlabels$itemnr))


# then remove them from dataframe, since we need a dataframe with only item data for the Rasch analyses
df <- df %>% select(all_of(c("Time",itemlabels$itemnr))) # add time here ? 
dfnotime <- df %>% select(!Time)
source("RISE_theme.R")

5 Demographics descriptives

Code
dif_spec <- dif %>% filter(Time=='SCREEN') %>% select(!all_of(c("Time","Patient")))
summary(dif_spec)
# RIdemographics(dif_spec,'Treatment') <- function crashes TODO make own 
 Treatment  sex           age        TreatmentAccessStart         education   
 MDD:2987   F:4005   Min.   :18.00   Min.   :2008         Primary      : 504  
 PD :1719   M:2456   1st Qu.:27.00   1st Qu.:2012         Secondary    :3152  
 SAD:1755            Median :33.00   Median :2014         Postsecondary:2805  
                     Mean   :35.38   Mean   :2014                             
                     3rd Qu.:42.00   3rd Qu.:2017                             
                     Max.   :84.00   Max.   :2019                             

6 Overall number of responses

Code
# Collapsed 
RIallresp(dfnotime)
# Seperate 
RIallresp_over_times <- df %>% # 
  split(.$Time) %>% # split the data 
  map(~ RIallresp(.x %>% dplyr::select(!Time))) #+ labs(title = .x$Time)) # create separate for each time

# make nice later 
RI_allresp_kable_grid = combine_kables_grid(RIallresp_over_times,cols=3)
RI_allresp_kable_grid
Response category
Number of responses
Percent
0
147379
24.9
1
148593
25.1
2
152831
25.8
3
67964
11.5
4
61598
10.4
5
10877
1.8
6
2301
0.4
x
SCREEN.Response category SCREEN.Number of responses SCREEN.Percent PRE.Response category PRE.Number of responses PRE.Percent WEEK01.Response category WEEK01.Number of responses WEEK01.Percent
0 8232 14.2 0 9700 17.3 0 9626 20.1
1 7983 13.7 1 9965 17.8 1 10868 22.7
2 15050 25.9 2 15616 27.9 2 13754 28.8
3 10126 17.4 3 9146 16.3 3 6674 14.0
4 13231 22.8 4 9629 17.2 4 5814 12.2
5 2816 4.8 5 1608 2.9 5 935 2.0
6 711 1.2 6 307 0.5 6 155 0.3
WEEK02.Response category WEEK02.Number of responses WEEK02.Percent WEEK03.Response category WEEK03.Number of responses WEEK03.Percent WEEK04.Response category WEEK04.Number of responses WEEK04.Percent
0 10426 21.2 0 10968 22.7 0 11626 24.7
1 11747 23.9 1 12365 25.6 1 12516 26.6
2 14076 28.6 2 13422 27.8 2 12690 26.9
3 6367 12.9 3 5755 11.9 3 5254 11.2
4 5502 11.2 4 4697 9.7 4 4155 8.8
5 884 1.8 5 856 1.8 5 722 1.5
6 183 0.4 6 150 0.3 6 134 0.3
WEEK05.Response category WEEK05.Number of responses WEEK05.Percent WEEK06.Response category WEEK06.Number of responses WEEK06.Percent WEEK07.Response category WEEK07.Number of responses WEEK07.Percent
0 11754 26.2 0 11745 27.3 0 11877 28.9
1 12299 27.4 1 12301 28.6 1 11965 29.1
2 11657 26.0 2 10933 25.4 2 10044 24.5
3 4732 10.6 3 4201 9.8 3 3801 9.3
4 3639 8.1 4 3242 7.5 4 2804 6.8
5 639 1.4 5 525 1.2 5 460 1.1
6 100 0.2 6 136 0.3 6 98 0.2
WEEK08.Response category WEEK08.Number of responses WEEK08.Percent WEEK09.Response category WEEK09.Number of responses WEEK09.Percent WEEK10.Response category WEEK10.Number of responses WEEK10.Percent
0 12015 30.6 0 11977 32.1 0 11951 34.0
1 11516 29.3 1 11197 30.0 1 10795 30.7
2 9632 24.5 2 8719 23.4 2 7755 22.0
3 3218 8.2 3 2913 7.8 3 2580 7.3
4 2431 6.2 4 2113 5.7 4 1775 5.0
5 408 1.0 5 335 0.9 5 272 0.8
6 92 0.2 6 69 0.2 6 71 0.2
POST.Response category POST.Number of responses POST.Percent " " " "
0 15482 34.9
1 13076 29.5
2 9483 21.4
3 3197 7.2
4 2566 5.8
5 417 0.9
6 95 0.2

7 Descriptives of raw data

Response distribution for all items are summarized below.

7.1 Floor/ceiling effects

Code
RIrawdist(dfnotime)

Code
#df_test = df %>% filter(Time=='SCREEN') %>% select(!Time)
# Across time 

# seperated time - might need to fix this remake in ggplot?
png(filename=paste0("./plots/",itemStart,"_rawdist_over_time.png"),
width=800,height=2400)
par(mfrow=c(7,2),mar = c(4, 4, 4, 2))  # c(5 <- bottom lines, 4, 4, 2)
RIrawdist_over_times <- df %>% # 
  split(.$Time) %>% # split the data 
  imap(~ {
    RIrawdist(.x %>% dplyr::select(!Time))
    mtext(paste("Time:", .y), side = 1, line = 3, adj = 0, cex = 0.8)
    })

dev.off()
png 
  2 
Code
print('Saved plot to /plots')
[1] "Saved plot to /plots"

7.2 Guttman structure

Code
RIheatmap(dfnotime)

Code
ggplots_with_spec_func_over_time(df,func_call = RIheatmap,title_all='Guttman structure over time',blankx=TRUE)

7.3 Descriptives - item level

Code
# over all times 
RItileplot(dfnotime)

Code
# by time 
ggplots_with_spec_func_over_time(df,func_call = RItileplot,title_all='',blankx=FALSE)

Code
RIbarstack(dfnotime)

Code
# by time 
ggplots_with_spec_func_over_time(df,func_call = RIbarstack,title_all='',blankx=FALSE)

Code
RIbarplot(dfnotime)

8 Rasch analysis 1

The eRm package, which uses Conditional Maximum Likelihood (CML) estimation, will be used primarily. For this analysis, the Partial Credit Model will be used.

Code
df_test = df %>% filter(Time=='SCREEN') %>% select(!Time)
RIlistItemsMargin(df, fontsize = 12)
itemnr item
MADRS.S_1 Sadness
MADRS.S_2 Anxiety
MADRS.S_3 Reduced sleep
MADRS.S_4 Reduced appetite
MADRS.S_5 Concentration difficulties
MADRS.S_6 Lassitude
MADRS.S_7 Inability to feel
MADRS.S_8 Pessimism
MADRS.S_9 Suicidal thoughts

8.1 Timepoint decision

Due to the amount of data across timepoints targeting will first be inspected, and the best fitting timepoint will be chosen based on targeting, which in turn will inform the rest of the analyses.

8.1.1 Targeting

Code
# increase fig-height above as needed, if you have many items og. value was 5. 
#RItargeting(dfnotime)

require(patchwork)
targeting_time_plots <- ggplots_with_spec_func_over_time(
  df,func_call = RItargeting,title_all='',
  blankx=FALSE,colsplot = 3,rowsplot=5,return_list_plots=TRUE)
good_layout_targeting_plots <- wrap_plots(targeting_time_plots, ncol = 3) + plot_layout(guides= "collect")
#print(good_layout_targeting_plots)
ggsave(
  filename = paste0("./plots/",items_to_use,"_targeting_over_time.png"),
  plot = good_layout_targeting_plots,
  width = 20,  # Increase width (in inches)
  height = 45, # Increase height (in inches)
  dpi = 300    # High resolution
)
print(good_layout_targeting_plots)

8.1.2 DIF treatment group targeting

Week02 looks promising.

Looking at treatment group for timepoints

Code
screenplot <- treatment_dif_by_time(timepoint='SCREEN',df=df,dif=dif,given_func=RItargeting)
wrap_plots(screenplot, ncol = 3)

Code
preplot <- treatment_dif_by_time(timepoint='PRE',df=df,dif=dif,given_func=RItargeting)
wrap_plots(preplot, ncol = 3)

Code
week01plot <- treatment_dif_by_time(timepoint='WEEK01',df=df,dif=dif,given_func=RItargeting)
wrap_plots(week01plot, ncol = 3)

Code
week02plot <- treatment_dif_by_time(timepoint='WEEK02',df=df,dif=dif,given_func=RItargeting)
wrap_plots(week02plot, ncol = 3)

8.2 Analysis continuation week02

Code
df_week02 <- df %>% filter(Time=='WEEK02') %>% select(!Time) #%>% sample_n(800)
Code
RIitemfit(df_week02, cutoff = "Smith98")
Item InfitMSQ Location 1 ± 2 / √n
MADRS.S_1 0.688 0.26 [0.973, 1.027]
MADRS.S_2 1.108 -0.78 [0.973, 1.027]
MADRS.S_3 1.593 -0.16 [0.973, 1.027]
MADRS.S_4 1.322 0.99 [0.973, 1.027]
MADRS.S_5 0.938 -0.17 [0.973, 1.027]
MADRS.S_6 0.793 -0.27 [0.973, 1.027]
MADRS.S_7 0.735 0.14 [0.973, 1.027]
MADRS.S_8 1.058 -0.53 [0.973, 1.027]
MADRS.S_9 0.846 0.60 [0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).
Code
simfit1 <- RIgetfit(df_week02, iterations = 300, cpu = n_cores)  

RIitemfit(df_week02, simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
MADRS.S_1 0.688 [0.95, 1.045] 0.677 [0.953, 1.054] 0.262 0.276 0.26
MADRS.S_2 1.108 [0.954, 1.046] 1.118 [0.953, 1.05] 0.062 0.068 -0.78
MADRS.S_3 1.593 [0.955, 1.052] 1.71 [0.939, 1.06] 0.541 0.650 -0.16
MADRS.S_4 1.322 [0.951, 1.056] 1.437 [0.923, 1.113] 0.266 0.324 0.99
MADRS.S_5 0.938 [0.955, 1.045] 0.929 [0.957, 1.051] 0.017 0.028 -0.17
MADRS.S_6 0.793 [0.954, 1.041] 0.788 [0.954, 1.042] 0.161 0.166 -0.27
MADRS.S_7 0.735 [0.955, 1.04] 0.725 [0.955, 1.04] 0.220 0.230 0.14
MADRS.S_8 1.058 [0.95, 1.052] 1.083 [0.946, 1.056] 0.006 0.027 -0.53
MADRS.S_9 0.846 [0.965, 1.049] 0.803 [0.955, 1.062] 0.119 0.152 0.60
Note:
MSQ values based on conditional calculations (n = 5465 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RIgetfitPlot(simfit1, df_week02)

Code
ICCplot(as.data.frame(df_week02), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[1:4])

[1] "Please press Zoom on the Plots window to see the plot"
Code
### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)
Code
ICCplot(as.data.frame(df_week02), 
        itemnumber = 5:8, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[5:8])

[1] "Please press Zoom on the Plots window to see the plot"
Code
RIrestscore(df_week02)
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
MADRS.S_1 0.75 0.61 0.14 0.000 *** 0.26 1.66
MADRS.S_2 0.59 0.62 0.03 0.000 *** -0.78 0.62
MADRS.S_3 0.44 0.63 0.19 0.000 *** -0.16 1.24
MADRS.S_4 0.46 0.61 0.15 0.000 *** 0.99 2.38
MADRS.S_5 0.65 0.62 0.03 0.000 *** -0.17 1.22
MADRS.S_6 0.72 0.62 0.10 0.000 *** -0.27 1.13
MADRS.S_7 0.74 0.60 0.14 0.000 *** 0.14 1.54
MADRS.S_8 0.63 0.64 0.01 0.055 . -0.53 0.87
MADRS.S_9 0.69 0.60 0.09 0.000 *** 0.60 1.99
Code
RIbootRestscore(df_week02,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600
Item Item-restscore result % of iterations Conditional MSQ infit Relative average item location
MADRS.S_1 overfit 100.0 0.69 1.66
MADRS.S_2 no misfit 82.4 1.11 0.62
MADRS.S_2 underfit 17.6 1.11 0.62
MADRS.S_3 underfit 100.0 1.59 1.24
MADRS.S_4 underfit 100.0 1.32 2.38
MADRS.S_5 no misfit 57.8 0.94 1.22
MADRS.S_5 overfit 42.2 0.94 1.22
MADRS.S_6 overfit 100.0 0.79 1.13
MADRS.S_7 overfit 100.0 0.74 1.54
MADRS.S_8 no misfit 95.0 1.06 0.87
MADRS.S_9 overfit 98.8 0.85 1.99
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.
Code
RIpcmPCA(df_week02)
PCA of Rasch model residuals
Eigenvalues Proportion of variance
1.70 22.5%
1.50 16.9%
1.31 14.5%
1.08 14.1%
0.97 10.4%
Code
simcor1 <- RIgetResidCor(df_week02, iterations = 500, cpu = n_cores)
RIresidcorr(df_week02, cutoff = simcor1$p99)
MADRS.S_1 MADRS.S_2 MADRS.S_3 MADRS.S_4 MADRS.S_5 MADRS.S_6 MADRS.S_7 MADRS.S_8 MADRS.S_9
MADRS.S_1
MADRS.S_2 -0.06
MADRS.S_3 -0.25 -0.12
MADRS.S_4 -0.16 -0.15 0
MADRS.S_5 -0.16 -0.08 -0.15 -0.13
MADRS.S_6 -0.05 -0.17 -0.26 -0.13 0.14
MADRS.S_7 0.15 -0.26 -0.24 -0.14 -0.06 0.09
MADRS.S_8 -0.09 -0.01 -0.26 -0.26 -0.2 -0.15 -0.09
MADRS.S_9 0.12 -0.21 -0.23 -0.12 -0.21 -0.11 0.09 0.07
Note:
Relative cut-off value is -0.046, which is 0.061 above the average correlation (-0.107).
Correlations above the cut-off are highlighted in red text.
Code
RIpartgamLD(df_week02)
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
MADRS.S_6 MADRS.S_5 0.349 0.019 0.310 0.387 0.000
MADRS.S_1 MADRS.S_7 0.338 0.021 0.296 0.380 0.000
MADRS.S_1 MADRS.S_9 0.334 0.021 0.292 0.375 0.000
MADRS.S_7 MADRS.S_1 0.325 0.021 0.283 0.367 0.000
MADRS.S_7 MADRS.S_9 0.315 0.022 0.271 0.358 0.000
MADRS.S_5 MADRS.S_6 0.305 0.020 0.266 0.344 0.000
MADRS.S_9 MADRS.S_1 0.294 0.022 0.252 0.337 0.000
MADRS.S_9 MADRS.S_7 0.294 0.022 0.250 0.337 0.000
MADRS.S_7 MADRS.S_6 0.291 0.021 0.249 0.333 0.000
MADRS.S_9 MADRS.S_8 0.277 0.020 0.239 0.316 0.000
MADRS.S_6 MADRS.S_7 0.272 0.022 0.230 0.314 0.000
MADRS.S_8 MADRS.S_9 0.214 0.020 0.174 0.253 0.000
MADRS.S_8 MADRS.S_2 0.172 0.019 0.135 0.209 0.000
MADRS.S_2 MADRS.S_8 0.162 0.018 0.126 0.198 0.000
MADRS.S_1 MADRS.S_8 0.13 0.020 0.089 0.170 0.000
MADRS.S_4 MADRS.S_3 0.119 0.020 0.081 0.158 0.000
MADRS.S_1 MADRS.S_2 0.11 0.021 0.068 0.152 0.000
MADRS.S_1 MADRS.S_6 0.104 0.022 0.061 0.147 0.000
MADRS.S_3 MADRS.S_4 0.104 0.019 0.066 0.143 0.000
MADRS.S_7 MADRS.S_5 0.095 0.022 0.052 0.139 0.001
MADRS.S_7 MADRS.S_8 0.085 0.021 0.043 0.126 0.004
MADRS.S_5 MADRS.S_3 0.072 0.019 0.035 0.109 0.009
MADRS.S_5 MADRS.S_2 0.07 0.020 0.031 0.110 0.031
Code
RIloadLoc(df_week02)

Code
mirt(df_week02, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code
# for fewer items or a more magnified figure, use:
#RIitemCats(df)
Code
# increase fig-height above as needed, if you have many items
RItargeting(df_week02)

Code
RIitemHierarchy(df_week02)

Code
iarm::score_groups(as.data.frame(df_week02)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)
  score_group    n
1           1 2752
2           2 2713
Code
dif_plots <- df_week02 %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code
clr_tests(df_week02, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr  df pvalue sig 
overall 1359 52 0       ***
Code
item_obsexp(PCM(df_week02))
Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.719    0.828   -8.015  -- 
MADRS.S_2  1.828    1.792    2.338  +  
MADRS.S_3  1.214    1.018   12.700  ++ 
MADRS.S_4  0.339    0.241   10.175  ++ 
MADRS.S_5  1.186    1.191   -0.344     
MADRS.S_6  1.111    1.180   -4.877  -- 
MADRS.S_7  0.912    1.006   -6.777  -- 
MADRS.S_8  1.426    1.416    0.608     
MADRS.S_9  0.577    0.639   -4.778  -- 

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1   2.430    2.323    5.836 ++ 
MADRS.S_2   3.293    3.329   -2.123 -  
MADRS.S_3   2.516    2.708  -10.544 -- 
MADRS.S_4   1.264    1.361   -5.422 -- 
MADRS.S_5   2.719    2.714    0.276    
MADRS.S_6   2.828    2.760    3.704 ++ 
MADRS.S_7   2.403    2.311    5.610 ++ 
MADRS.S_8   3.369    3.379   -0.559    
MADRS.S_9   1.897    1.835    4.007 ++ 
Code
grouping_based_on_score = score_groups(as.data.frame(df_week02))
partgam_DIF(dat.items = as.data.frame(df_week02),
            dat.exo = grouping_based_on_score) 
       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_2 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_3 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 MADRS.S_4 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 MADRS.S_5 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
6 MADRS.S_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
7 MADRS.S_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
8 MADRS.S_8 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
9 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
Code
dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_week02, dif.sex)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.260 0.243 0.251 0.013 0.018
MADRS.S_2 -0.834 -0.723 -0.779 0.078 0.111
MADRS.S_3 -0.135 -0.254 -0.194 0.084 0.119
MADRS.S_4 0.907 1.155 1.031 0.175 0.248
MADRS.S_5 -0.155 -0.258 -0.207 0.073 0.103
MADRS.S_6 -0.196 -0.433 -0.314 0.168 0.237
MADRS.S_7 0.192 0.004 0.098 0.133 0.188
MADRS.S_8 -0.594 -0.439 -0.516 0.110 0.155
MADRS.S_9 0.554 0.705 0.629 0.106 0.150
Code
partgam_DIF(dat.items = as.data.frame(df_week02),
            dat.exo = dif.sex) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.sex  0.2073 0.0263 0.0000  0.0000  ***  0.1558  0.2588
2 MADRS.S_2 dif.sex -0.1617 0.0239 0.0000  0.0000  *** -0.2085 -0.1148
3 MADRS.S_3 dif.sex  0.0010 0.0228 0.9641  1.0000      -0.0438  0.0458
4 MADRS.S_4 dif.sex -0.1020 0.0274 0.0002  0.0018   ** -0.1558 -0.0482
5 MADRS.S_5 dif.sex  0.0548 0.0254 0.0312  0.2811       0.0049  0.1046
6 MADRS.S_6 dif.sex  0.0786 0.0260 0.0025  0.0225    *  0.0277  0.1296
7 MADRS.S_7 dif.sex  0.2151 0.0266 0.0000  0.0000  ***  0.1630  0.2672
8 MADRS.S_8 dif.sex -0.1738 0.0233 0.0000  0.0000  *** -0.2195 -0.1280
9 MADRS.S_9 dif.sex  0.0961 0.0276 0.0005  0.0045   **  0.0420  0.1502
Code
dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_week02, dif.age)

Item 3 4 6 8 9 Mean location StDev MaxDiff
MADRS.S_1 0.201 0.220 0.308 0.272 -0.077 0.185 0.152 0.385
MADRS.S_2 -0.873 -0.405 -0.902 -0.444 -0.773 -0.679 0.238 0.497
MADRS.S_3 0.039 0.050 -0.279 -0.356 -0.192 -0.148 0.185 0.406
MADRS.S_4 0.854 0.544 1.000 0.902 0.687 0.798 0.181 0.456
MADRS.S_5 -0.160 -0.210 -0.229 -0.204 0.135 -0.133 0.152 0.364
MADRS.S_6 -0.315 -0.240 -0.274 -0.285 -0.430 -0.309 0.073 0.190
MADRS.S_7 0.149 0.118 0.217 -0.003 0.262 0.149 0.102 0.265
MADRS.S_8 -0.484 -0.305 -0.564 -0.476 -0.167 -0.399 0.161 0.397
MADRS.S_9 0.589 0.228 0.722 0.595 0.555 0.538 0.185 0.494
Code
partgam_DIF(dat.items = as.data.frame(df_week02),
            dat.exo = dif.age) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.0380 0.0158 0.0161  0.1452      -0.0690 -0.0070
2 MADRS.S_2 dif.age -0.1415 0.0140 0.0000  0.0000  *** -0.1689 -0.1142
3 MADRS.S_3 dif.age  0.1160 0.0132 0.0000  0.0000  ***  0.0902  0.1418
4 MADRS.S_4 dif.age -0.0367 0.0158 0.0201  0.1808      -0.0676 -0.0058
5 MADRS.S_5 dif.age  0.0134 0.0149 0.3690  1.0000      -0.0158  0.0427
6 MADRS.S_6 dif.age -0.0317 0.0152 0.0369  0.3317      -0.0614 -0.0019
7 MADRS.S_7 dif.age  0.1092 0.0158 0.0000  0.0000  ***  0.0782  0.1401
8 MADRS.S_8 dif.age -0.0270 0.0137 0.0485  0.4366      -0.0539 -0.0002
9 MADRS.S_9 dif.age  0.0690 0.0161 0.0000  0.0002  ***  0.0374  0.1006
Code
dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_week02, dif.edu)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.271 0.210 0.241 0.043 0.061
MADRS.S_2 -0.792 -0.782 -0.787 0.007 0.010
MADRS.S_3 -0.159 -0.178 -0.169 0.014 0.020
MADRS.S_4 0.882 1.169 1.026 0.203 0.287
MADRS.S_5 -0.206 -0.154 -0.180 0.037 0.053
MADRS.S_6 -0.333 -0.190 -0.261 0.101 0.143
MADRS.S_7 0.144 0.109 0.126 0.025 0.035
MADRS.S_8 -0.454 -0.678 -0.566 0.158 0.224
MADRS.S_9 0.647 0.493 0.570 0.109 0.154
Code
partgam_DIF(dat.items = as.data.frame(df_week02),
            dat.exo = dif.edu) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.edu  0.0241 0.0238 0.3127  1.0000      -0.0226  0.0708
2 MADRS.S_2 dif.edu -0.0422 0.0215 0.0495  0.4453      -0.0844 -0.0001
3 MADRS.S_3 dif.edu  0.0231 0.0202 0.2529  1.0000      -0.0165  0.0627
4 MADRS.S_4 dif.edu -0.0631 0.0239 0.0082  0.0734    . -0.1099 -0.0164
5 MADRS.S_5 dif.edu -0.0262 0.0228 0.2487  1.0000      -0.0708  0.0183
6 MADRS.S_6 dif.edu -0.0888 0.0232 0.0001  0.0012   ** -0.1343 -0.0432
7 MADRS.S_7 dif.edu  0.0704 0.0243 0.0037  0.0333    *  0.0229  0.1180
8 MADRS.S_8 dif.edu  0.0871 0.0214 0.0000  0.0004  ***  0.0451  0.1290
9 MADRS.S_9 dif.edu  0.0530 0.0247 0.0320  0.2878       0.0046  0.1015

8.3 DIF-Year

Code
dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_week02, dif.yr)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.329 0.213 0.271 0.082 0.116
MADRS.S_2 -0.778 -0.807 -0.792 0.021 0.029
MADRS.S_3 -0.172 -0.175 -0.174 0.002 0.003
MADRS.S_4 1.051 0.946 0.998 0.074 0.105
MADRS.S_5 -0.088 -0.240 -0.164 0.107 0.152
MADRS.S_6 -0.319 -0.253 -0.286 0.047 0.066
MADRS.S_7 0.066 0.182 0.124 0.083 0.117
MADRS.S_8 -0.557 -0.534 -0.546 0.016 0.023
MADRS.S_9 0.469 0.669 0.569 0.141 0.200
Code
partgam_DIF(dat.items = as.data.frame(df_week02),
            dat.exo = dif.yr) 
       Item    Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.yr -0.0132 0.0165 0.4233  1.0000      -0.0455  0.0191
2 MADRS.S_2 dif.yr  0.0585 0.0151 0.0001  0.0010   **  0.0288  0.0881
3 MADRS.S_3 dif.yr -0.0031 0.0142 0.8297  1.0000      -0.0310  0.0249
4 MADRS.S_4 dif.yr -0.0271 0.0169 0.1084  0.9757      -0.0602  0.0060
5 MADRS.S_5 dif.yr  0.0094 0.0157 0.5488  1.0000      -0.0213  0.0402
6 MADRS.S_6 dif.yr -0.0044 0.0160 0.7812  1.0000      -0.0358  0.0269
7 MADRS.S_7 dif.yr -0.0177 0.0168 0.2921  1.0000      -0.0507  0.0152
8 MADRS.S_8 dif.yr  0.0228 0.0147 0.1222  1.0000      -0.0061  0.0517
9 MADRS.S_9 dif.yr -0.0704 0.0168 0.0000  0.0002  *** -0.1033 -0.0375

8.4 Analysis part 1 decision

Remove item 3 to begin with : motivation goes here

9 Rasch analysis 2

Code
df_r <- df_week02 %>% select(!MADRS.S_3) 
Code
RIlistItemsMargin(df_r, fontsize = 12)
itemnr item
MADRS.S_1 Sadness
MADRS.S_2 Anxiety
MADRS.S_4 Reduced appetite
MADRS.S_5 Concentration difficulties
MADRS.S_6 Lassitude
MADRS.S_7 Inability to feel
MADRS.S_8 Pessimism
MADRS.S_9 Suicidal thoughts
Code
RIitemfit(df_r, cutoff = "Smith98")
Item InfitMSQ Location 1 ± 2 / √n
MADRS.S_1 0.719 0.27 [0.973, 1.027]
MADRS.S_2 1.234 -0.88 [0.973, 1.027]
MADRS.S_4 1.501 1.06 [0.973, 1.027]
MADRS.S_5 1.03 -0.20 [0.973, 1.027]
MADRS.S_6 0.831 -0.30 [0.973, 1.027]
MADRS.S_7 0.771 0.14 [0.973, 1.027]
MADRS.S_8 1.131 -0.60 [0.973, 1.027]
MADRS.S_9 0.891 0.61 [0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).
Code
simfit1 <- RIgetfit(df_r, iterations = 300, cpu = n_cores)  

RIitemfit(df_r, simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
MADRS.S_1 0.719 [0.944, 1.048] 0.71 [0.94, 1.046] 0.225 0.23 0.27
MADRS.S_2 1.234 [0.954, 1.046] 1.254 [0.954, 1.049] 0.188 0.205 -0.88
MADRS.S_4 1.501 [0.943, 1.059] 1.811 [0.906, 1.179] 0.442 0.632 1.06
MADRS.S_5 1.03 [0.96, 1.039] 1.023 [0.959, 1.038] no misfit no misfit -0.20
MADRS.S_6 0.831 [0.952, 1.048] 0.832 [0.953, 1.047] 0.121 0.121 -0.30
MADRS.S_7 0.771 [0.956, 1.046] 0.758 [0.953, 1.051] 0.185 0.195 0.14
MADRS.S_8 1.131 [0.951, 1.054] 1.16 [0.95, 1.056] 0.077 0.104 -0.60
MADRS.S_9 0.891 [0.951, 1.04] 0.842 [0.942, 1.054] 0.06 0.1 0.61
Note:
MSQ values based on conditional calculations (n = 5465 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RIgetfitPlot(simfit1, df_r)

Code
ICCplot(as.data.frame(df_r), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[1:4])

[1] "Please press Zoom on the Plots window to see the plot"
Code
### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)
Code
ICCplot(as.data.frame(df_r), 
        itemnumber = 5:8, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[5:8])

[1] "Please press Zoom on the Plots window to see the plot"
Code
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
MADRS.S_1 0.76 0.65 0.11 0.000 *** 0.27 1.79
MADRS.S_2 0.59 0.66 0.07 0.000 *** -0.88 0.65
MADRS.S_4 0.45 0.64 0.19 0.000 *** 1.06 2.59
MADRS.S_5 0.65 0.65 0.00 0.844 -0.20 1.32
MADRS.S_6 0.73 0.65 0.08 0.000 *** -0.30 1.22
MADRS.S_7 0.75 0.64 0.11 0.000 *** 0.14 1.67
MADRS.S_8 0.64 0.67 0.03 0.000 *** -0.60 0.93
MADRS.S_9 0.70 0.64 0.06 0.000 *** 0.61 2.14
Code
RIbootRestscore(df_r,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600
Item Item-restscore result % of iterations Conditional MSQ infit Relative average item location
MADRS.S_1 overfit 100.0 0.72 1.79
MADRS.S_2 no misfit 10.2 1.23 0.65
MADRS.S_2 underfit 89.8 1.23 0.65
MADRS.S_4 underfit 100.0 1.50 2.59
MADRS.S_5 no misfit 98.0 1.03 1.32
MADRS.S_6 overfit 99.4 0.83 1.22
MADRS.S_7 overfit 100.0 0.77 1.67
MADRS.S_8 no misfit 67.8 1.13 0.93
MADRS.S_8 underfit 32.2 1.13 0.93
MADRS.S_9 no misfit 10.6 0.89 2.14
MADRS.S_9 overfit 89.4 0.89 2.14
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.
Code
RIpcmPCA(df_r)
PCA of Rasch model residuals
Eigenvalues Proportion of variance
1.58 23.8%
1.53 18.8%
1.19 17.8%
1.11 12.9%
0.90 9.9%
Code
simcor1 <- RIgetResidCor(df_r, iterations = 500, cpu = n_cores)
RIresidcorr(df_r, cutoff = simcor1$p99)
MADRS.S_1 MADRS.S_2 MADRS.S_4 MADRS.S_5 MADRS.S_6 MADRS.S_7 MADRS.S_8 MADRS.S_9
MADRS.S_1
MADRS.S_2 -0.1
MADRS.S_4 -0.18 -0.15
MADRS.S_5 -0.21 -0.11 -0.13
MADRS.S_6 -0.13 -0.22 -0.15 0.1
MADRS.S_7 0.09 -0.31 -0.16 -0.11 0.03
MADRS.S_8 -0.16 -0.04 -0.28 -0.26 -0.23 -0.17
MADRS.S_9 0.06 -0.25 -0.13 -0.25 -0.18 0.04 0.02
Note:
Relative cut-off value is -0.066, which is 0.061 above the average correlation (-0.127).
Correlations above the cut-off are highlighted in red text.
Code
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
MADRS.S_6 MADRS.S_5 0.343 0.020 0.305 0.382 0.000
MADRS.S_1 MADRS.S_7 0.295 0.022 0.251 0.339 0.000
MADRS.S_1 MADRS.S_9 0.293 0.022 0.249 0.336 0.000
MADRS.S_5 MADRS.S_6 0.283 0.021 0.243 0.323 0.000
MADRS.S_7 MADRS.S_1 0.283 0.022 0.240 0.327 0.000
MADRS.S_7 MADRS.S_9 0.278 0.023 0.233 0.322 0.000
MADRS.S_7 MADRS.S_6 0.249 0.022 0.206 0.293 0.000
MADRS.S_9 MADRS.S_7 0.247 0.023 0.201 0.292 0.000
MADRS.S_9 MADRS.S_1 0.245 0.022 0.201 0.288 0.000
MADRS.S_9 MADRS.S_8 0.239 0.020 0.199 0.280 0.000
MADRS.S_6 MADRS.S_7 0.226 0.023 0.182 0.270 0.000
MADRS.S_8 MADRS.S_9 0.168 0.021 0.127 0.208 0.000
MADRS.S_8 MADRS.S_2 0.162 0.019 0.126 0.199 0.000
MADRS.S_2 MADRS.S_8 0.137 0.019 0.100 0.173 0.000
MADRS.S_1 MADRS.S_2 0.107 0.022 0.064 0.150 0.000
MADRS.S_1 MADRS.S_4 0.096 0.023 0.051 0.141 0.002
MADRS.S_6 MADRS.S_4 0.092 0.022 0.049 0.136 0.002
MADRS.S_7 MADRS.S_4 0.087 0.023 0.041 0.133 0.011
MADRS.S_7 MADRS.S_5 0.083 0.023 0.038 0.127 0.015
MADRS.S_5 MADRS.S_2 0.08 0.020 0.041 0.119 0.003
MADRS.S_1 MADRS.S_8 0.074 0.021 0.033 0.116 0.025
Code
RIloadLoc(df_r)

Code
mirt(df_r, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code
# for fewer items or a more magnified figure, use:
#RIitemCats(df)
Code
# increase fig-height above as needed, if you have many items
RItargeting(df_r)

Code

Code
iarm::score_groups(as.data.frame(df_r)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)
  score_group    n
1           1 2709
2           2 2756
Code
dif_plots <- df_r %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code
clr_tests(df_r, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr df pvalue sig 
overall 967 46 0       ***
Code
Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.709    0.792   -6.435  -- 
MADRS.S_2  1.826    1.758    4.655  ++ 
MADRS.S_4  0.347    0.209   15.423  ++ 
MADRS.S_5  1.172    1.157    1.107     
MADRS.S_6  1.098    1.145   -3.484  -- 
MADRS.S_7  0.892    0.972   -6.087  -- 
MADRS.S_8  1.401    1.371    1.921     
MADRS.S_9  0.559    0.600   -3.222  -- 

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1  2.418    2.340    4.542  ++ 
MADRS.S_2  3.286    3.349   -4.018  -- 
MADRS.S_4  1.244    1.375   -7.717  -- 
MADRS.S_5  2.717    2.731   -0.854     
MADRS.S_6  2.821    2.776    2.573  +  
MADRS.S_7  2.406    2.330    4.886  ++ 
MADRS.S_8  3.373    3.402   -1.700     
MADRS.S_9  1.897    1.858    2.620  +  
Code
grouping_based_on_score = score_groups(as.data.frame(df_r))
partgam_DIF(dat.items = as.data.frame(df_r),
            dat.exo = grouping_based_on_score) 
       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_2 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_4 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 MADRS.S_5 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 MADRS.S_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
6 MADRS.S_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
7 MADRS.S_8 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
8 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
Code
dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_r, dif.sex)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.270 0.226 0.248 0.031 0.044
MADRS.S_2 -0.934 -0.832 -0.883 0.072 0.102
MADRS.S_4 0.974 1.240 1.107 0.188 0.266
MADRS.S_5 -0.183 -0.308 -0.246 0.089 0.125
MADRS.S_6 -0.225 -0.493 -0.359 0.189 0.268
MADRS.S_7 0.196 -0.015 0.090 0.149 0.211
MADRS.S_8 -0.668 -0.521 -0.594 0.104 0.147
MADRS.S_9 0.571 0.704 0.637 0.094 0.133
Code
partgam_DIF(dat.items = as.data.frame(df_r),
            dat.exo = dif.sex) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.sex  0.2209 0.0265 0.0000  0.0000  ***  0.1689  0.2728
2 MADRS.S_2 dif.sex -0.1625 0.0239 0.0000  0.0000  *** -0.2093 -0.1156
3 MADRS.S_4 dif.sex -0.1092 0.0273 0.0001  0.0005  *** -0.1626 -0.0557
4 MADRS.S_5 dif.sex  0.0543 0.0255 0.0330  0.2641       0.0044  0.1042
5 MADRS.S_6 dif.sex  0.0838 0.0263 0.0014  0.0116    *  0.0322  0.1354
6 MADRS.S_7 dif.sex  0.2262 0.0269 0.0000  0.0000  ***  0.1735  0.2788
7 MADRS.S_8 dif.sex -0.1769 0.0236 0.0000  0.0000  *** -0.2232 -0.1306
8 MADRS.S_9 dif.sex  0.1004 0.0280 0.0003  0.0026   **  0.0456  0.1552
Code
dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_r, dif.age)

Item 3 4 6 8 9 Mean location StDev MaxDiff
MADRS.S_1 0.281 0.114 -0.216 0.265 -0.001 0.088 0.206 0.497
MADRS.S_2 -0.966 -1.008 -0.862 -0.672 -0.773 -0.856 0.138 0.336
MADRS.S_4 0.883 1.040 1.075 1.004 0.745 0.950 0.135 0.330
MADRS.S_5 -0.228 -0.028 -0.205 -0.248 0.228 -0.096 0.201 0.476
MADRS.S_6 -0.332 -0.330 -0.203 -0.332 -0.373 -0.314 0.065 0.170
MADRS.S_7 0.123 0.382 0.193 -0.010 -0.396 0.058 0.291 0.778
MADRS.S_8 -0.510 -0.740 -0.474 -0.617 -0.161 -0.500 0.216 0.578
MADRS.S_9 0.749 0.570 0.692 0.611 0.731 0.671 0.078 0.180
Code
partgam_DIF(dat.items = as.data.frame(df_r),
            dat.exo = dif.age) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.0067 0.0162 0.6794  1.0000      -0.0384  0.0250
2 MADRS.S_2 dif.age -0.1251 0.0140 0.0000  0.0000  *** -0.1526 -0.0976
3 MADRS.S_4 dif.age -0.0255 0.0156 0.1038  0.8303      -0.0561  0.0052
4 MADRS.S_5 dif.age  0.0386 0.0149 0.0095  0.0760    .  0.0094  0.0677
5 MADRS.S_6 dif.age -0.0020 0.0154 0.8975  1.0000      -0.0321  0.0281
6 MADRS.S_7 dif.age  0.1424 0.0159 0.0000  0.0000  ***  0.1113  0.1736
7 MADRS.S_8 dif.age -0.0070 0.0139 0.6141  1.0000      -0.0342  0.0202
8 MADRS.S_9 dif.age  0.0941 0.0163 0.0000  0.0000  ***  0.0622  0.1260
Code
dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_r, dif.edu)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.272 0.213 0.242 0.041 0.058
MADRS.S_2 -0.890 -0.887 -0.889 0.002 0.003
MADRS.S_4 0.951 1.241 1.096 0.205 0.290
MADRS.S_5 -0.241 -0.187 -0.214 0.038 0.054
MADRS.S_6 -0.374 -0.222 -0.298 0.107 0.152
MADRS.S_7 0.144 0.102 0.123 0.030 0.042
MADRS.S_8 -0.522 -0.761 -0.641 0.169 0.239
MADRS.S_9 0.661 0.502 0.581 0.112 0.159
Code
partgam_DIF(dat.items = as.data.frame(df_r),
            dat.exo = dif.edu) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.edu  0.0287 0.0243 0.2381  1.0000      -0.0190  0.0763
2 MADRS.S_2 dif.edu -0.0383 0.0214 0.0738  0.5902      -0.0802  0.0037
3 MADRS.S_4 dif.edu -0.0618 0.0238 0.0095  0.0761    . -0.1085 -0.0151
4 MADRS.S_5 dif.edu -0.0232 0.0229 0.3096  1.0000      -0.0681  0.0216
5 MADRS.S_6 dif.edu -0.0912 0.0235 0.0001  0.0008  *** -0.1371 -0.0452
6 MADRS.S_7 dif.edu  0.0730 0.0247 0.0031  0.0246    *  0.0247  0.1214
7 MADRS.S_8 dif.edu  0.0894 0.0216 0.0000  0.0003  ***  0.0471  0.1317
8 MADRS.S_9 dif.edu  0.0582 0.0250 0.0198  0.1580       0.0093  0.1071
Code
dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_r, dif.yr)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_r),
            dat.exo = dif.yr) 
       Item    Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.yr -0.0159 0.0167 0.3418  1.0000      -0.0486  0.0168
2 MADRS.S_2 dif.yr  0.0562 0.0150 0.0002  0.0015   **  0.0267  0.0856
3 MADRS.S_4 dif.yr -0.0241 0.0168 0.1511  1.0000      -0.0569  0.0088
4 MADRS.S_5 dif.yr  0.0086 0.0157 0.5856  1.0000      -0.0222  0.0394
5 MADRS.S_6 dif.yr -0.0024 0.0163 0.8849  1.0000      -0.0342  0.0295
6 MADRS.S_7 dif.yr -0.0183 0.0170 0.2818  1.0000      -0.0517  0.0151
7 MADRS.S_8 dif.yr  0.0219 0.0149 0.1399  1.0000      -0.0072  0.0510
8 MADRS.S_9 dif.yr -0.0754 0.0170 0.0000  0.0001  *** -0.1087 -0.0421

9.1 Analysis part 2 decision

Remove item 4, appetite : from the residual contrast loadings we see it loads to something completely different than the others.

10 Rasch analysis 3

Removing item 4

Code
df_r3 <- df_week02 %>% select(!c(MADRS.S_4,MADRS.S_3))
Code
RIlistItemsMargin(df_r3, fontsize = 12)
itemnr item
MADRS.S_1 Sadness
MADRS.S_2 Anxiety
MADRS.S_5 Concentration difficulties
MADRS.S_6 Lassitude
MADRS.S_7 Inability to feel
MADRS.S_8 Pessimism
MADRS.S_9 Suicidal thoughts
Code
RIitemfit(df_r3, cutoff = "Smith98")
Item InfitMSQ Location 1 ± 2 / √n
MADRS.S_1 0.775 0.45 [0.973, 1.027]
MADRS.S_2 1.325 -0.77 [0.973, 1.027]
MADRS.S_5 1.119 -0.05 [0.973, 1.027]
MADRS.S_6 0.901 -0.15 [0.973, 1.027]
MADRS.S_7 0.823 0.34 [0.973, 1.027]
MADRS.S_8 1.169 -0.49 [0.973, 1.027]
MADRS.S_9 0.958 0.81 [0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).
Code
simfit1 <- RIgetfit(df_r3, iterations = 300, cpu = n_cores)  

RIitemfit(df_r3, simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
MADRS.S_1 0.775 [0.952, 1.049] 0.76 [0.949, 1.06] 0.177 0.189 0.45
MADRS.S_2 1.325 [0.963, 1.038] 1.342 [0.959, 1.041] 0.287 0.301 -0.77
MADRS.S_5 1.119 [0.955, 1.047] 1.105 [0.956, 1.046] 0.072 0.059 -0.05
MADRS.S_6 0.901 [0.952, 1.052] 0.896 [0.954, 1.046] 0.051 0.058 -0.15
MADRS.S_7 0.823 [0.954, 1.043] 0.807 [0.952, 1.046] 0.131 0.145 0.34
MADRS.S_8 1.169 [0.95, 1.046] 1.194 [0.95, 1.048] 0.123 0.146 -0.49
MADRS.S_9 0.958 [0.955, 1.051] 0.908 [0.938, 1.073] no misfit 0.030 0.81
Note:
MSQ values based on conditional calculations (n = 5465 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RIgetfitPlot(simfit1, df_r3)

Code
ICCplot(as.data.frame(df_r3), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[1:4])

[1] "Please press Zoom on the Plots window to see the plot"
Code
### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)
Code
ICCplot(as.data.frame(df_r3), 
        itemnumber = 5:7, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[5:7])

[1] "Please press Zoom on the Plots window to see the plot"
Code
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
MADRS.S_1 0.76 0.67 0.09 0.000 *** 0.45 1.94
MADRS.S_2 0.59 0.68 0.09 0.000 *** -0.77 0.71
MADRS.S_5 0.65 0.68 0.03 0.002 ** -0.05 1.44
MADRS.S_6 0.72 0.68 0.04 0.000 *** -0.15 1.33
MADRS.S_7 0.75 0.67 0.08 0.000 *** 0.34 1.82
MADRS.S_8 0.65 0.70 0.05 0.000 *** -0.49 0.99
MADRS.S_9 0.70 0.66 0.04 0.000 *** 0.81 2.29
Code
RIbootRestscore(df_r3,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600
Item Item-restscore result % of iterations Conditional MSQ infit Relative average item location
MADRS.S_1 overfit 99.8 0.78 1.94
MADRS.S_2 underfit 99.8 1.33 0.71
MADRS.S_5 no misfit 90.0 1.12 1.44
MADRS.S_5 underfit 10.0 1.12 1.44
MADRS.S_6 no misfit 25.4 0.90 1.33
MADRS.S_6 overfit 74.6 0.90 1.33
MADRS.S_7 overfit 99.8 0.82 1.82
MADRS.S_8 no misfit 52.0 1.17 0.99
MADRS.S_8 underfit 48.0 1.17 0.99
MADRS.S_9 no misfit 62.0 0.96 2.29
MADRS.S_9 overfit 38.0 0.96 2.29
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.
Code
RIpcmPCA(df_r3)
PCA of Rasch model residuals
Eigenvalues Proportion of variance
1.61 24.7%
1.57 23.5%
1.18 17%
0.93 12.9%
0.88 11.4%
Code
simcor1 <- RIgetResidCor(df_r3, iterations = 500, cpu = n_cores)
RIresidcorr(df_r3, cutoff = simcor1$p99)
MADRS.S_1 MADRS.S_2 MADRS.S_5 MADRS.S_6 MADRS.S_7 MADRS.S_8 MADRS.S_9
MADRS.S_1
MADRS.S_2 -0.13
MADRS.S_5 -0.25 -0.13
MADRS.S_6 -0.17 -0.25 0.08
MADRS.S_7 0.06 -0.34 -0.14 0
MADRS.S_8 -0.22 -0.08 -0.31 -0.28 -0.22
MADRS.S_9 0.03 -0.28 -0.28 -0.21 0.01 -0.02
Note:
Relative cut-off value is -0.086, which is 0.063 above the average correlation (-0.149).
Correlations above the cut-off are highlighted in red text.
Code
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
MADRS.S_6 MADRS.S_5 0.346 0.020 0.308 0.385 0.000
MADRS.S_1 MADRS.S_7 0.298 0.022 0.254 0.341 0.000
MADRS.S_1 MADRS.S_9 0.295 0.022 0.252 0.338 0.000
MADRS.S_7 MADRS.S_1 0.284 0.022 0.240 0.327 0.000
MADRS.S_5 MADRS.S_6 0.281 0.020 0.241 0.321 0.000
MADRS.S_7 MADRS.S_9 0.275 0.023 0.230 0.319 0.000
MADRS.S_7 MADRS.S_6 0.25 0.022 0.207 0.293 0.000
MADRS.S_9 MADRS.S_1 0.243 0.022 0.199 0.286 0.000
MADRS.S_9 MADRS.S_7 0.242 0.023 0.197 0.287 0.000
MADRS.S_6 MADRS.S_7 0.23 0.022 0.186 0.274 0.000
MADRS.S_9 MADRS.S_8 0.208 0.021 0.167 0.248 0.000
MADRS.S_8 MADRS.S_9 0.141 0.021 0.100 0.183 0.000
MADRS.S_8 MADRS.S_2 0.138 0.019 0.101 0.176 0.000
MADRS.S_2 MADRS.S_8 0.101 0.019 0.064 0.139 0.000
MADRS.S_1 MADRS.S_2 0.099 0.022 0.056 0.142 0.000
MADRS.S_7 MADRS.S_5 0.084 0.023 0.039 0.128 0.009
MADRS.S_5 MADRS.S_2 0.07 0.020 0.031 0.109 0.021
Code
RIloadLoc(df_r3)

Code
mirt(df_r3, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code
# for fewer items or a more magnified figure, use:
#RIitemCats(df)
Code
# increase fig-height above as needed, if you have many items
RItargeting(df_r3)

Code

Code
iarm::score_groups(as.data.frame(df_r3)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)
  score_group    n
1           1 2898
2           2 2567
Code
dif_plots <- df_r3 %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code
clr_tests(df_r3, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr df pvalue sig 
overall 595 40 0       ***
Code
Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.755    0.817   -5.036  -- 
MADRS.S_2  1.866    1.792    5.391  ++ 
MADRS.S_5  1.223    1.187    2.773  +  
MADRS.S_6  1.152    1.173   -1.635     
MADRS.S_7  0.945    0.999   -4.325  -- 
MADRS.S_8  1.459    1.418    2.722  +  
MADRS.S_9  0.608    0.621   -1.092     

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1  2.495    2.428    3.940  ++ 
MADRS.S_2  3.353    3.433   -5.121  -- 
MADRS.S_5  2.777    2.816   -2.356  -  
MADRS.S_6  2.891    2.868    1.343     
MADRS.S_7  2.462    2.403    3.823  ++ 
MADRS.S_8  3.458    3.502   -2.705  -  
MADRS.S_9  1.943    1.929    0.989     
Code
grouping_based_on_score = score_groups(as.data.frame(df_r3))
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = grouping_based_on_score) 
       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_2 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_5 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 MADRS.S_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 MADRS.S_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
6 MADRS.S_8 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
7 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
Code
dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_r3, dif.sex)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.443 0.423 0.433 0.014 0.020
MADRS.S_2 -0.855 -0.693 -0.774 0.114 0.162
MADRS.S_5 -0.044 -0.129 -0.087 0.060 0.085
MADRS.S_6 -0.089 -0.325 -0.207 0.166 0.235
MADRS.S_7 0.375 0.191 0.283 0.130 0.184
MADRS.S_8 -0.584 -0.377 -0.480 0.146 0.207
MADRS.S_9 0.754 0.909 0.831 0.110 0.155
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.sex) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.sex  0.1993 0.0267 0.0000  0.0000  ***  0.1469  0.2517
2 MADRS.S_2 dif.sex -0.1754 0.0239 0.0000  0.0000  *** -0.2222 -0.1286
3 MADRS.S_5 dif.sex  0.0365 0.0255 0.1526  1.0000      -0.0135  0.0865
4 MADRS.S_6 dif.sex  0.0685 0.0264 0.0095  0.0667    .  0.0167  0.1202
5 MADRS.S_7 dif.sex  0.2089 0.0271 0.0000  0.0000  ***  0.1558  0.2619
6 MADRS.S_8 dif.sex -0.2011 0.0237 0.0000  0.0000  *** -0.2476 -0.1547
7 MADRS.S_9 dif.sex  0.0903 0.0280 0.0013  0.0088   **  0.0354  0.1452
Code
dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_r3, dif.age)

Item 3 4 6 8 9 Mean location StDev MaxDiff
MADRS.S_1 0.428 0.301 -0.098 0.449 0.016 0.219 0.248 0.547
MADRS.S_2 -0.885 -0.915 -0.784 -0.479 -0.885 -0.789 0.181 0.436
MADRS.S_5 -0.101 0.120 -0.018 -0.139 0.255 0.023 0.163 0.394
MADRS.S_6 -0.208 -0.183 -0.026 -0.233 -0.428 -0.216 0.144 0.403
MADRS.S_7 0.277 0.563 0.408 0.146 0.384 0.356 0.156 0.417
MADRS.S_8 -0.414 -0.648 -0.405 -0.537 -0.102 -0.421 0.204 0.546
MADRS.S_9 0.905 0.760 0.923 0.794 0.761 0.828 0.079 0.163
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.age) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.0134 0.0162 0.4054  1.0000      -0.0451  0.0182
2 MADRS.S_2 dif.age -0.1291 0.0140 0.0000  0.0000  *** -0.1565 -0.1016
3 MADRS.S_5 dif.age  0.0339 0.0148 0.0222  0.1555       0.0049  0.0630
4 MADRS.S_6 dif.age -0.0081 0.0154 0.5984  1.0000      -0.0384  0.0221
5 MADRS.S_7 dif.age  0.1370 0.0159 0.0000  0.0000  ***  0.1058  0.1682
6 MADRS.S_8 dif.age -0.0097 0.0140 0.4882  1.0000      -0.0371  0.0177
7 MADRS.S_9 dif.age  0.0904 0.0163 0.0000  0.0000  ***  0.0584  0.1223
Code
dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_r3, dif.edu)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.433 0.427 0.430 0.004 0.006
MADRS.S_2 -0.802 -0.768 -0.785 0.024 0.034
MADRS.S_5 -0.107 -0.008 -0.057 0.070 0.099
MADRS.S_6 -0.248 -0.043 -0.146 0.145 0.205
MADRS.S_7 0.312 0.320 0.316 0.006 0.008
MADRS.S_8 -0.424 -0.648 -0.536 0.159 0.225
MADRS.S_9 0.835 0.720 0.778 0.082 0.115
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.edu) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.edu  0.0136 0.0243 0.5752  1.0000      -0.0340  0.0613
2 MADRS.S_2 dif.edu -0.0427 0.0214 0.0464  0.3247      -0.0848 -0.0007
3 MADRS.S_5 dif.edu -0.0330 0.0228 0.1478  1.0000      -0.0777  0.0117
4 MADRS.S_6 dif.edu -0.1015 0.0234 0.0000  0.0001  *** -0.1474 -0.0556
5 MADRS.S_7 dif.edu  0.0671 0.0246 0.0063  0.0439    *  0.0190  0.1153
6 MADRS.S_8 dif.edu  0.0909 0.0218 0.0000  0.0002  ***  0.0482  0.1336
7 MADRS.S_9 dif.edu  0.0493 0.0249 0.0477  0.3338       0.0005  0.0982
Code
dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_r3, dif.yr)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.yr) 
       Item    Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.yr -0.0195 0.0167 0.2413  1.0000      -0.0522  0.0131
2 MADRS.S_2 dif.yr  0.0509 0.0151 0.0007  0.0052   **  0.0213  0.0805
3 MADRS.S_5 dif.yr  0.0069 0.0157 0.6613  1.0000      -0.0239  0.0376
4 MADRS.S_6 dif.yr -0.0063 0.0163 0.7005  1.0000      -0.0383  0.0257
5 MADRS.S_7 dif.yr -0.0182 0.0170 0.2862  1.0000      -0.0515  0.0152
6 MADRS.S_8 dif.yr  0.0190 0.0151 0.2063  1.0000      -0.0105  0.0486
7 MADRS.S_9 dif.yr -0.0780 0.0170 0.0000  0.0000  *** -0.1113 -0.0448

10.1 Analysis part 3 decision

High number of boot misfit for item 2 and for conditional item fit it was most deviant - removing.

11 Rasch analysis 4

Removing item 2

Code
df_r3 <- df_week02 %>% select(!c(MADRS.S_4,MADRS.S_3,MADRS.S_2))
Code
RIlistItemsMargin(df_r3, fontsize = 12)
itemnr item
MADRS.S_1 Sadness
MADRS.S_5 Concentration difficulties
MADRS.S_6 Lassitude
MADRS.S_7 Inability to feel
MADRS.S_8 Pessimism
MADRS.S_9 Suicidal thoughts
Code
RIitemfit(df_r3, cutoff = "Smith98")
Item InfitMSQ Location 1 ± 2 / √n
MADRS.S_1 0.838 0.34 [0.973, 1.027]
MADRS.S_5 1.215 -0.19 [0.973, 1.027]
MADRS.S_6 0.929 -0.30 [0.973, 1.027]
MADRS.S_7 0.806 0.22 [0.973, 1.027]
MADRS.S_8 1.321 -0.66 [0.973, 1.027]
MADRS.S_9 0.967 0.70 [0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).
Code
simfit1 <- RIgetfit(df_r3, iterations = 300, cpu = n_cores)  

RIitemfit(df_r3, simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
MADRS.S_1 0.838 [0.95, 1.05] 0.843 [0.943, 1.054] 0.112 0.100 0.34
MADRS.S_5 1.215 [0.949, 1.059] 1.214 [0.947, 1.057] 0.156 0.157 -0.19
MADRS.S_6 0.929 [0.953, 1.04] 0.93 [0.954, 1.038] 0.024 0.024 -0.30
MADRS.S_7 0.806 [0.955, 1.048] 0.793 [0.955, 1.05] 0.149 0.162 0.22
MADRS.S_8 1.321 [0.955, 1.057] 1.358 [0.954, 1.057] 0.264 0.301 -0.66
MADRS.S_9 0.967 [0.959, 1.049] 0.923 [0.951, 1.067] no misfit 0.028 0.70
Note:
MSQ values based on conditional calculations (n = 5465 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RIgetfitPlot(simfit1, df_r3)

Code
ICCplot(as.data.frame(df_r3), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[1:4])

[1] "Please press Zoom on the Plots window to see the plot"
Code
### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)
Code
ICCplot(as.data.frame(df_r3), 
        itemnumber = 5:6, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[5:6])

[1] "Please press Zoom on the Plots window to see the plot"
Code
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
MADRS.S_1 0.76 0.69 0.07 0 *** 0.34 2.03
MADRS.S_5 0.64 0.70 0.06 0 *** -0.19 1.50
MADRS.S_6 0.73 0.70 0.03 0 *** -0.30 1.40
MADRS.S_7 0.77 0.69 0.08 0 *** 0.22 1.91
MADRS.S_8 0.64 0.71 0.07 0 *** -0.66 1.03
MADRS.S_9 0.72 0.69 0.03 0 *** 0.70 2.40
Code
RIbootRestscore(df_r3,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600
Item Item-restscore result % of iterations Conditional MSQ infit Relative average item location
MADRS.S_1 overfit 96.0 0.84 2.03
MADRS.S_5 no misfit 25.6 1.21 1.50
MADRS.S_5 underfit 74.4 1.21 1.50
MADRS.S_6 no misfit 53.4 0.93 1.40
MADRS.S_6 overfit 46.6 0.93 1.40
MADRS.S_7 overfit 100.0 0.81 1.91
MADRS.S_8 underfit 98.6 1.32 1.03
MADRS.S_9 no misfit 68.8 0.97 2.40
MADRS.S_9 overfit 31.2 0.97 2.40
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.
Code
RIpcmPCA(df_r3)
PCA of Rasch model residuals
Eigenvalues Proportion of variance
1.66 30.8%
1.37 23.5%
1.03 16.7%
0.97 15.1%
0.95 13.7%
Code
simcor1 <- RIgetResidCor(df_r3, iterations = 500, cpu = n_cores)
RIresidcorr(df_r3, cutoff = simcor1$p99)
MADRS.S_1 MADRS.S_5 MADRS.S_6 MADRS.S_7 MADRS.S_8 MADRS.S_9
MADRS.S_1
MADRS.S_5 -0.27
MADRS.S_6 -0.22 0.04
MADRS.S_7 0 -0.2 -0.09
MADRS.S_8 -0.23 -0.32 -0.32 -0.28
MADRS.S_9 -0.02 -0.34 -0.29 -0.07 -0.06
Note:
Relative cut-off value is -0.114, which is 0.064 above the average correlation (-0.178).
Correlations above the cut-off are highlighted in red text.
Code
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
MADRS.S_6 MADRS.S_5 0.363 0.020 0.324 0.401 0
MADRS.S_5 MADRS.S_6 0.266 0.021 0.226 0.306 0
MADRS.S_7 MADRS.S_1 0.258 0.023 0.213 0.303 0
MADRS.S_1 MADRS.S_9 0.247 0.023 0.202 0.292 0
MADRS.S_1 MADRS.S_7 0.239 0.023 0.194 0.285 0
MADRS.S_9 MADRS.S_8 0.222 0.021 0.180 0.263 0
MADRS.S_9 MADRS.S_1 0.212 0.023 0.167 0.257 0
MADRS.S_7 MADRS.S_6 0.19 0.023 0.144 0.236 0
MADRS.S_7 MADRS.S_9 0.186 0.025 0.138 0.235 0
MADRS.S_9 MADRS.S_7 0.135 0.025 0.087 0.184 0
MADRS.S_6 MADRS.S_7 0.133 0.024 0.087 0.180 0
MADRS.S_8 MADRS.S_9 0.1 0.022 0.057 0.142 0
Code
RIloadLoc(df_r3)

Code
mirt(df_r3, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code
# for fewer items or a more magnified figure, use:
#RIitemCats(df)
Code
# increase fig-height above as needed, if you have many items
RItargeting(df_r3)

Code

Code
iarm::score_groups(as.data.frame(df_r3)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)
  score_group    n
1           1 2798
2           2 2667
Code
dif_plots <- df_r3 %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code
clr_tests(df_r3, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr df pvalue sig 
overall 523 34 0       ***
Code
Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.744    0.799   -4.514  -- 
MADRS.S_5  1.234    1.181    4.225  ++ 
MADRS.S_6  1.150    1.166   -1.252     
MADRS.S_7  0.937    0.988   -4.135  -- 
MADRS.S_8  1.478    1.399    5.501  ++ 
MADRS.S_9  0.583    0.595   -0.968     

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1  2.466    2.413    3.336  ++ 
MADRS.S_5  2.750    2.802   -3.360  -- 
MADRS.S_6  2.867    2.852    0.952     
MADRS.S_7  2.445    2.395    3.441  ++ 
MADRS.S_8  3.414    3.492   -5.020  -- 
MADRS.S_9  1.938    1.927    0.826     
Code
grouping_based_on_score = score_groups(as.data.frame(df_r3))
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = grouping_based_on_score) 
       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_5 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 MADRS.S_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 MADRS.S_8 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
6 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
Code
dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_r3, dif.sex)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.323 0.324 0.324 0.001 0.001
MADRS.S_5 -0.196 -0.255 -0.225 0.042 0.059
MADRS.S_6 -0.243 -0.457 -0.350 0.152 0.214
MADRS.S_7 0.246 0.086 0.166 0.114 0.161
MADRS.S_8 -0.770 -0.520 -0.645 0.177 0.251
MADRS.S_9 0.639 0.821 0.730 0.129 0.182
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.sex) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.sex  0.1585 0.0272 0.0000  0.0000  ***  0.1052  0.2117
2 MADRS.S_5 dif.sex  0.0048 0.0256 0.8518  1.0000      -0.0454  0.0549
3 MADRS.S_6 dif.sex  0.0207 0.0270 0.4420  1.0000      -0.0321  0.0736
4 MADRS.S_7 dif.sex  0.1807 0.0282 0.0000  0.0000  ***  0.1254  0.2360
5 MADRS.S_8 dif.sex -0.2368 0.0234 0.0000  0.0000  *** -0.2826 -0.1909
6 MADRS.S_9 dif.sex  0.0619 0.0287 0.0310  0.1859       0.0057  0.1181
Code
dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_r3, dif.age)

Item 3 4 6 7 Mean location StDev MaxDiff
MADRS.S_1 0.289 0.161 0.456 -0.188 0.180 0.273 0.644
MADRS.S_5 -0.251 -0.033 -0.273 -0.006 -0.141 0.141 0.267
MADRS.S_6 -0.362 -0.348 -0.306 -0.045 -0.265 0.149 0.317
MADRS.S_7 0.137 0.427 0.100 0.104 0.192 0.157 0.327
MADRS.S_8 -0.581 -0.839 -0.703 -0.442 -0.641 0.170 0.397
MADRS.S_9 0.769 0.632 0.726 0.577 0.676 0.087 0.192
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.age) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.0446 0.0162 0.0059  0.0355    * -0.0764 -0.0128
2 MADRS.S_5 dif.age  0.0128 0.0148 0.3903  1.0000      -0.0163  0.0419
3 MADRS.S_6 dif.age -0.0379 0.0156 0.0153  0.0920    . -0.0685 -0.0073
4 MADRS.S_7 dif.age  0.1152 0.0165 0.0000  0.0000  ***  0.0828  0.1475
5 MADRS.S_8 dif.age -0.0370 0.0140 0.0082  0.0493    * -0.0645 -0.0096
6 MADRS.S_9 dif.age  0.0624 0.0167 0.0002  0.0012   **  0.0296  0.0952
Code
dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_r3, dif.edu)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.318 0.323 0.321 0.003 0.005
MADRS.S_5 -0.251 -0.144 -0.198 0.076 0.108
MADRS.S_6 -0.397 -0.181 -0.289 0.153 0.216
MADRS.S_7 0.191 0.207 0.199 0.011 0.016
MADRS.S_8 -0.590 -0.827 -0.708 0.167 0.237
MADRS.S_9 0.728 0.621 0.675 0.076 0.108
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.edu) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.edu -0.0010 0.0246 0.9691  1.0000      -0.0492  0.0473
2 MADRS.S_5 dif.edu -0.0392 0.0228 0.0864  0.5182      -0.0839  0.0056
3 MADRS.S_6 dif.edu -0.1139 0.0238 0.0000  0.0000  *** -0.1606 -0.0672
4 MADRS.S_7 dif.edu  0.0589 0.0256 0.0213  0.1279       0.0088  0.1091
5 MADRS.S_8 dif.edu  0.0770 0.0218 0.0004  0.0025   **  0.0342  0.1198
6 MADRS.S_9 dif.edu  0.0442 0.0256 0.0842  0.5050      -0.0060  0.0944
Code
dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_r3, dif.yr)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.yr) 
       Item    Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.yr -0.0037 0.0169 0.8269  1.0000      -0.0367  0.0294
2 MADRS.S_5 dif.yr  0.0177 0.0157 0.2589  1.0000      -0.0130  0.0485
3 MADRS.S_6 dif.yr  0.0030 0.0166 0.8575  1.0000      -0.0296  0.0356
4 MADRS.S_7 dif.yr -0.0105 0.0177 0.5536  1.0000      -0.0452  0.0242
5 MADRS.S_8 dif.yr  0.0297 0.0150 0.0480  0.2880       0.0003  0.0591
6 MADRS.S_9 dif.yr -0.0679 0.0173 0.0001  0.0005  *** -0.1019 -0.0340

11.1 Analysis part 4 decision

Due to the general misfit of infit and outfit (which was the largest) for item 8, aswell as its item-rest-score bootstrap underfitting, removing this.

12 Rasch analysis 5

Removing item 8

Code
df_r3 <- df %>% filter(Time=='WEEK02') %>% select(!Time) %>% select(!c(MADRS.S_4,MADRS.S_3,MADRS.S_2,MADRS.S_8))
Code
RIlistItemsMargin(df_r3, fontsize = 12)
itemnr item
MADRS.S_1 Sadness
MADRS.S_5 Concentration difficulties
MADRS.S_6 Lassitude
MADRS.S_7 Inability to feel
MADRS.S_9 Suicidal thoughts
Code
RIitemfit(df_r3, cutoff = "Smith98")
Item InfitMSQ Location 1 ± 2 / √n
MADRS.S_1 0.889 0.22 [0.973, 1.027]
MADRS.S_5 1.253 -0.35 [0.973, 1.027]
MADRS.S_6 0.943 -0.46 [0.973, 1.027]
MADRS.S_7 0.825 0.08 [0.973, 1.027]
MADRS.S_9 1.113 0.59 [0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).
Code
simfit1 <- RIgetfit(df_r3, iterations = 300, cpu = n_cores)  

RIitemfit(df_r3, simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
MADRS.S_1 0.889 [0.953, 1.051] 0.912 [0.952, 1.046] 0.064 0.04 0.22
MADRS.S_5 1.253 [0.96, 1.046] 1.254 [0.962, 1.041] 0.207 0.213 -0.35
MADRS.S_6 0.943 [0.943, 1.044] 0.947 [0.944, 1.043] no misfit no misfit -0.46
MADRS.S_7 0.825 [0.959, 1.048] 0.82 [0.957, 1.046] 0.134 0.137 0.08
MADRS.S_9 1.113 [0.955, 1.039] 1.093 [0.94, 1.057] 0.074 0.036 0.59
Note:
MSQ values based on conditional calculations (n = 5465 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RIgetfitPlot(simfit1, df_r3)

Code
ICCplot(as.data.frame(df_r3), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[c(1,5,6,7)])

[1] "Please press Zoom on the Plots window to see the plot"
Code
### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)
Code
ICCplot(as.data.frame(df_r3), 
        itemnumber = c(5), 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[9])

Code
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
MADRS.S_1 0.75 0.72 0.03 0.000 *** 0.22 2.14
MADRS.S_5 0.66 0.72 0.06 0.000 *** -0.35 1.57
MADRS.S_6 0.75 0.72 0.03 0.000 *** -0.46 1.46
MADRS.S_7 0.78 0.71 0.07 0.000 *** 0.08 2.00
MADRS.S_9 0.69 0.71 0.02 0.039 * 0.59 2.51
Code
RIbootRestscore(df_r3,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600
Item Item-restscore result % of iterations Conditional MSQ infit Relative average item location
MADRS.S_1 no misfit 41.8 0.89 2.14
MADRS.S_1 overfit 58.2 0.89 2.14
MADRS.S_5 no misfit 12.2 1.25 1.57
MADRS.S_5 underfit 87.8 1.25 1.57
MADRS.S_6 no misfit 73.2 0.94 1.46
MADRS.S_6 overfit 26.8 0.94 1.46
MADRS.S_7 overfit 98.4 0.82 2.00
MADRS.S_9 no misfit 96.4 1.11 2.51
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.
Code
RIpcmPCA(df_r3)
PCA of Rasch model residuals
Eigenvalues Proportion of variance
1.72 36.8%
1.16 23.2%
1.07 21.1%
1.05 18.9%
0.01 0.1%
Code
simcor1 <- RIgetResidCor(df_r3, iterations = 500, cpu = n_cores)
RIresidcorr(df_r3, cutoff = simcor1$p99)
MADRS.S_1 MADRS.S_5 MADRS.S_6 MADRS.S_7 MADRS.S_9
MADRS.S_1
MADRS.S_5 -0.37
MADRS.S_6 -0.32 -0.06
MADRS.S_7 -0.08 -0.32 -0.2
MADRS.S_9 -0.03 -0.39 -0.34 -0.11
Note:
Relative cut-off value is -0.163, which is 0.059 above the average correlation (-0.222).
Correlations above the cut-off are highlighted in red text.
Code
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
MADRS.S_6 MADRS.S_5 0.324 0.021 0.282 0.366 0.000
MADRS.S_1 MADRS.S_9 0.312 0.023 0.267 0.356 0.000
MADRS.S_7 MADRS.S_9 0.243 0.025 0.194 0.292 0.000
MADRS.S_9 MADRS.S_1 0.242 0.023 0.198 0.287 0.000
MADRS.S_7 MADRS.S_1 0.229 0.024 0.181 0.276 0.000
MADRS.S_5 MADRS.S_6 0.188 0.022 0.145 0.232 0.000
MADRS.S_1 MADRS.S_7 0.176 0.025 0.127 0.225 0.000
MADRS.S_9 MADRS.S_7 0.125 0.025 0.076 0.175 0.000
MADRS.S_7 MADRS.S_6 0.087 0.026 0.037 0.138 0.014
Code
RIloadLoc(df_r3)

Code
mirt(df_r3, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code
# for fewer items or a more magnified figure, use:
#RIitemCats(df)
Code
# increase fig-height above as needed, if you have many items
RItargeting(df_r3)

Code

Code
iarm::score_groups(as.data.frame(df_r3)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)
  score_group    n
1           1 2941
2           2 2524
Code
dif_plots <- df_r3 %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code
clr_tests(df_r3, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr df pvalue sig 
overall 460 28 0       ***
Code
Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.801    0.836   -3.024  -- 
MADRS.S_5  1.293    1.235    4.918  ++ 
MADRS.S_6  1.211    1.216   -0.407     
MADRS.S_7  1.001    1.033   -2.780  -  
MADRS.S_9  0.635    0.622    1.190     

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1  2.532    2.496    2.286  +  
MADRS.S_5  2.821    2.882   -3.951  -- 
MADRS.S_6  2.944    2.939    0.316     
MADRS.S_7  2.498    2.465    2.348  +  
MADRS.S_9  1.982    1.996   -1.058     
Code
grouping_based_on_score = score_groups(as.data.frame(df_r3))
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = grouping_based_on_score) 
       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_5 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 MADRS.S_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
Code
dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_r3, dif.sex)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.190 0.232 0.211 0.030 0.042
MADRS.S_5 -0.375 -0.371 -0.373 0.003 0.004
MADRS.S_6 -0.424 -0.579 -0.501 0.110 0.155
MADRS.S_7 0.098 -0.018 0.040 0.082 0.116
MADRS.S_9 0.511 0.736 0.624 0.159 0.224
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.sex) 
       Item     Var   gamma     se pvalue padj.BH sig   lower   upper
1 MADRS.S_1 dif.sex  0.0873 0.0280 0.0019  0.0093  **  0.0323  0.1422
2 MADRS.S_5 dif.sex -0.0571 0.0262 0.0296  0.1479     -0.1085 -0.0057
3 MADRS.S_6 dif.sex -0.0655 0.0281 0.0195  0.0976   . -0.1206 -0.0105
4 MADRS.S_7 dif.sex  0.0987 0.0298 0.0009  0.0046  **  0.0404  0.1571
5 MADRS.S_9 dif.sex -0.0121 0.0287 0.6726  1.0000     -0.0683  0.0441
Code
dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_r3, dif.age)

Item 3 4 5 Mean location StDev MaxDiff
MADRS.S_1 0.208 0.108 0.283 0.200 0.088 0.175
MADRS.S_5 -0.342 -0.296 -0.428 -0.356 0.067 0.132
MADRS.S_6 -0.421 -0.522 -0.470 -0.471 0.050 0.101
MADRS.S_7 0.017 0.145 0.029 0.064 0.070 0.128
MADRS.S_9 0.539 0.566 0.585 0.563 0.023 0.047
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.age) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.0565 0.0164 0.0006  0.0029   ** -0.0886 -0.0243
2 MADRS.S_5 dif.age -0.0001 0.0153 0.9966  1.0000      -0.0301  0.0299
3 MADRS.S_6 dif.age -0.0556 0.0163 0.0007  0.0034   ** -0.0876 -0.0235
4 MADRS.S_7 dif.age  0.1119 0.0173 0.0000  0.0000  ***  0.0780  0.1458
5 MADRS.S_9 dif.age  0.0491 0.0168 0.0034  0.0172    *  0.0162  0.0819
Code
dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_r3, dif.edu)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.218 0.176 0.197 0.029 0.041
MADRS.S_5 -0.391 -0.330 -0.360 0.042 0.060
MADRS.S_6 -0.542 -0.365 -0.453 0.125 0.177
MADRS.S_7 0.078 0.042 0.060 0.025 0.036
MADRS.S_9 0.637 0.478 0.557 0.113 0.160
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.edu) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.edu  0.0279 0.0251 0.2671  1.0000      -0.0214  0.0772
2 MADRS.S_5 dif.edu -0.0233 0.0235 0.3210  1.0000      -0.0692  0.0227
3 MADRS.S_6 dif.edu -0.1050 0.0248 0.0000  0.0001  *** -0.1536 -0.0564
4 MADRS.S_7 dif.edu  0.0928 0.0265 0.0005  0.0023   **  0.0408  0.1448
5 MADRS.S_9 dif.edu  0.0589 0.0257 0.0219  0.1096       0.0085  0.1092
Code
dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_r3, dif.yr)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.yr) 
       Item    Var   gamma     se pvalue padj.BH sig   lower   upper
1 MADRS.S_1 dif.yr  0.0144 0.0173 0.4043  1.0000     -0.0195  0.0484
2 MADRS.S_5 dif.yr  0.0318 0.0162 0.0501  0.2504      0.0000  0.0636
3 MADRS.S_6 dif.yr  0.0165 0.0173 0.3402  1.0000     -0.0174  0.0505
4 MADRS.S_7 dif.yr  0.0008 0.0183 0.9653  1.0000     -0.0351  0.0367
5 MADRS.S_9 dif.yr -0.0607 0.0174 0.0005  0.0024  ** -0.0947 -0.0267

12.1 Analysis part 5 decision

5 and 6 seems to load to something else and 5 has the worst fit (infit outfit). Also the partial gamma in relation to 5 and 6. While 7 underfits in the boostrap, the infit/outfit is not as deviant as 5. Removing 5.

13 Rasch analysis 6

Removing item 5

Code
df_r3 <- df %>% filter(Time=='WEEK02') %>% select(!Time) %>% 
  select(!c(MADRS.S_4,MADRS.S_3,MADRS.S_2,MADRS.S_8,MADRS.S_5))
Code
RIlistItemsMargin(df_r3, fontsize = 12)
itemnr item
MADRS.S_1 Sadness
MADRS.S_6 Lassitude
MADRS.S_7 Inability to feel
MADRS.S_9 Suicidal thoughts
Code
RIitemfit(df_r3, cutoff = "Smith98")
Item InfitMSQ Location 1 ± 2 / √n
MADRS.S_1 0.892 0.15 [0.973, 1.027]
MADRS.S_6 1.184 -0.57 [0.973, 1.027]
MADRS.S_7 0.843 -0.01 [0.973, 1.027]
MADRS.S_9 1.11 0.52 [0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).
Code
simfit1 <- RIgetfit(df_r3, iterations = 300, cpu = n_cores)  

RIitemfit(df_r3, simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
MADRS.S_1 0.892 [0.953, 1.043] 0.912 [0.95, 1.044] 0.061 0.038 0.15
MADRS.S_6 1.184 [0.959, 1.042] 1.196 [0.962, 1.044] 0.142 0.152 -0.57
MADRS.S_7 0.843 [0.953, 1.057] 0.837 [0.952, 1.054] 0.110 0.115 -0.01
MADRS.S_9 1.11 [0.96, 1.047] 1.096 [0.949, 1.062] 0.063 0.034 0.52
Note:
MSQ values based on conditional calculations (n = 5465 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RIgetfitPlot(simfit1, df_r3)

Code
ICCplot(as.data.frame(df_r3), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[c(1,6,7,9)])

[1] "Please press Zoom on the Plots window to see the plot"
Code
### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)
Code
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
MADRS.S_1 0.77 0.74 0.03 0.000 *** 0.15 2.27
MADRS.S_6 0.71 0.75 0.04 0.000 *** -0.57 1.55
MADRS.S_7 0.80 0.74 0.06 0.000 *** -0.01 2.11
MADRS.S_9 0.72 0.74 0.02 0.042 * 0.52 2.64
Code
RIbootRestscore(df_r3,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600
Item Item-restscore result % of iterations Conditional MSQ infit Relative average item location
MADRS.S_1 no misfit 59.8 0.89 2.27
MADRS.S_1 overfit 40.2 0.89 2.27
MADRS.S_6 no misfit 75.0 1.18 1.55
MADRS.S_6 underfit 25.0 1.18 1.55
MADRS.S_7 no misfit 7.8 0.84 2.11
MADRS.S_7 overfit 92.2 0.84 2.11
MADRS.S_9 no misfit 98.0 1.11 2.64
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.
Code
RIpcmPCA(df_r3)
PCA of Rasch model residuals
Eigenvalues Proportion of variance
1.51 41.7%
1.26 30.7%
1.22 27.4%
0.01 0.2%
Code
simcor1 <- RIgetResidCor(df_r3, iterations = 500, cpu = n_cores)
RIresidcorr(df_r3, cutoff = simcor1$p99)
MADRS.S_1 MADRS.S_6 MADRS.S_7 MADRS.S_9
MADRS.S_1
MADRS.S_6 -0.39
MADRS.S_7 -0.23 -0.25
MADRS.S_9 -0.19 -0.43 -0.26
Note:
Relative cut-off value is -0.225, which is 0.065 above the average correlation (-0.29).
Correlations above the cut-off are highlighted in red text.
Code
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
MADRS.S_7 MADRS.S_6 0.247 0.025 0.198 0.297 0.000
MADRS.S_1 MADRS.S_9 0.216 0.025 0.167 0.265 0.000
MADRS.S_7 MADRS.S_1 0.129 0.027 0.076 0.182 0.000
MADRS.S_7 MADRS.S_9 0.117 0.028 0.063 0.172 0.000
MADRS.S_9 MADRS.S_1 0.092 0.026 0.041 0.142 0.004
Code
RIloadLoc(df_r3)

Code
mirt(df_r3, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code
# for fewer items or a more magnified figure, use:
#RIitemCats(df)
Code
# increase fig-height above as needed, if you have many items
RItargeting(df_r3)

Code

Code
iarm::score_groups(as.data.frame(df_r3)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)
  score_group    n
1           1 2965
2           2 2500
Code
dif_plots <- df_r3 %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code
clr_tests(df_r3, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr df pvalue sig 
overall 316 22 0       ***
Code
Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.812    0.842   -2.778  -  
MADRS.S_6  1.270    1.235    3.255  ++ 
MADRS.S_7  1.026    1.047   -1.968     
MADRS.S_9  0.634    0.618    1.500     

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1  2.57     2.54     2.12   +  
MADRS.S_6  2.94     2.98    -2.50   -  
MADRS.S_7  2.52     2.50     1.65      
MADRS.S_9  2.02     2.04    -1.33      
Code
grouping_based_on_score = score_groups(as.data.frame(df_r3))
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = grouping_based_on_score) 
       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
Code
dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_r3, dif.sex)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.111 0.149 0.130 0.027 0.038
MADRS.S_6 -0.550 -0.686 -0.618 0.097 0.136
MADRS.S_7 0.003 -0.112 -0.055 0.082 0.116
MADRS.S_9 0.436 0.650 0.543 0.152 0.214
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.sex) 
       Item     Var   gamma     se pvalue padj.BH sig   lower   upper
1 MADRS.S_1 dif.sex  0.0734 0.0294 0.0125  0.0500   .  0.0158  0.1310
2 MADRS.S_6 dif.sex -0.0837 0.0278 0.0026  0.0105   * -0.1383 -0.0292
3 MADRS.S_7 dif.sex  0.0884 0.0311 0.0044  0.0177   *  0.0275  0.1493
4 MADRS.S_9 dif.sex -0.0391 0.0298 0.1891  0.7564     -0.0974  0.0193
Code
dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_r3, dif.age)

Item 3 4 5 Mean location StDev MaxDiff
MADRS.S_1 0.128 0.049 0.190 0.122 0.071 0.142
MADRS.S_6 -0.527 -0.632 -0.597 -0.585 0.054 0.105
MADRS.S_7 -0.063 0.071 -0.084 -0.025 0.084 0.155
MADRS.S_9 0.462 0.513 0.490 0.489 0.025 0.051
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.age) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.0688 0.0172 0.0001  0.0003  *** -0.1025 -0.0351
2 MADRS.S_6 dif.age -0.0539 0.0162 0.0009  0.0035   ** -0.0856 -0.0222
3 MADRS.S_7 dif.age  0.1161 0.0179 0.0000  0.0000  ***  0.0810  0.1511
4 MADRS.S_9 dif.age  0.0522 0.0175 0.0029  0.0117    *  0.0178  0.0866
Code
dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_r3, dif.edu)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.133 0.109 0.121 0.017 0.023
MADRS.S_6 -0.670 -0.469 -0.570 0.142 0.201
MADRS.S_7 -0.019 -0.050 -0.034 0.022 0.032
MADRS.S_9 0.556 0.410 0.483 0.103 0.146
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.edu) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.edu  0.0217 0.0263 0.4104  1.0000      -0.0299  0.0732
2 MADRS.S_6 dif.edu -0.1160 0.0245 0.0000  0.0000  *** -0.1640 -0.0679
3 MADRS.S_7 dif.edu  0.0811 0.0276 0.0033  0.0132    *  0.0270  0.1351
4 MADRS.S_9 dif.edu  0.0713 0.0267 0.0076  0.0304    *  0.0190  0.1237
Code
dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_r3, dif.yr)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.yr) 
       Item    Var   gamma     se pvalue padj.BH sig   lower   upper
1 MADRS.S_1 dif.yr  0.0235 0.0179 0.1899  0.7595     -0.0116  0.0587
2 MADRS.S_6 dif.yr  0.0293 0.0172 0.0887  0.3548     -0.0044  0.0631
3 MADRS.S_7 dif.yr  0.0063 0.0189 0.7403  1.0000     -0.0308  0.0434
4 MADRS.S_9 dif.yr -0.0611 0.0181 0.0007  0.0029  ** -0.0966 -0.0257

13.1 Analysis part 6 decision

7 or 6, trying to remove both in separate analyses.

14 Rasch analysis 7 part 1 of 2

Removing item 6

Code
df_r3 <- df %>% filter(Time=='WEEK02') %>% select(!Time) %>% 
  select(!c(MADRS.S_4,MADRS.S_3,MADRS.S_2,MADRS.S_8,MADRS.S_5,MADRS.S_6))
Code
RIlistItemsMargin(df_r3, fontsize = 12)
itemnr item
MADRS.S_1 Sadness
MADRS.S_7 Inability to feel
MADRS.S_9 Suicidal thoughts
Code
RIitemfit(df_r3, cutoff = "Smith98")
Item InfitMSQ Location 1 ± 2 / √n
MADRS.S_1 0.943 -0.05 [0.973, 1.027]
MADRS.S_7 0.971 -0.22 [0.973, 1.027]
MADRS.S_9 1.096 0.33 [0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).
Code
simfit1 <- RIgetfit(df_r3, iterations = 300, cpu = n_cores)  

RIitemfit(df_r3, simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
MADRS.S_1 0.943 [0.96, 1.039] 0.956 [0.96, 1.04] 0.017 0.004 -0.05
MADRS.S_7 0.971 [0.967, 1.033] 0.967 [0.966, 1.032] no misfit no misfit -0.22
MADRS.S_9 1.096 [0.956, 1.033] 1.086 [0.954, 1.038] 0.063 0.048 0.33
Note:
MSQ values based on conditional calculations (n = 5465 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RIgetfitPlot(simfit1, df_r3)

Code
ICCplot(as.data.frame(df_r3), 
        itemnumber = 1:3, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[c(1,7,9)])

[1] "Please press Zoom on the Plots window to see the plot"
Code
### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)
Code
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
MADRS.S_1 0.78 0.76 0.02 0.023 * -0.05 2.39
MADRS.S_7 0.78 0.76 0.02 0.044 * -0.22 2.22
MADRS.S_9 0.75 0.76 0.01 0.119 0.33 2.77
Code
RIbootRestscore(df_r3,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600
Item Item-restscore result % of iterations Conditional MSQ infit Relative average item location
MADRS.S_1 no misfit 96.2 0.94 2.39
MADRS.S_7 no misfit 98.4 0.97 2.22
MADRS.S_9 no misfit 100.0 1.10 2.77
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.
Code
RIpcmPCA(df_r3)
PCA of Rasch model residuals
Eigenvalues Proportion of variance
1.51 53.1%
1.49 46.6%
0.01 0.3%
Code
simcor1 <- RIgetResidCor(df_r3, iterations = 500, cpu = n_cores)
RIresidcorr(df_r3, cutoff = simcor1$p99)
MADRS.S_1 MADRS.S_7 MADRS.S_9
MADRS.S_1
MADRS.S_7 -0.37
MADRS.S_9 -0.43 -0.43
Note:
Relative cut-off value is -0.354, which is 0.06 above the average correlation (-0.413).
Correlations above the cut-off are highlighted in red text.
Code
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
MADRS.S_1 MADRS.S_7 0.11 0.03 0.051 0.170 0.002
MADRS.S_7 MADRS.S_1 0.103 0.03 0.044 0.162 0.004
Code
RIloadLoc(df_r3)

Code
mirt(df_r3, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code
# for fewer items or a more magnified figure, use:
#RIitemCats(df)
Code
# increase fig-height above as needed, if you have many items
RItargeting(df_r3)

Code

Code
iarm::score_groups(as.data.frame(df_r3)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)
  score_group    n
1           1 2961
2           2 2504
Code
dif_plots <- df_r3 %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code
clr_tests(df_r3, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr df pvalue sig 
overall 134 16 0       ***
Code
Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.868    0.894   -2.438  -  
MADRS.S_7  1.129    1.124    0.507     
MADRS.S_9  0.661    0.641    1.924     

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1  2.592    2.568    1.927     
MADRS.S_7  2.526    2.531   -0.423     
MADRS.S_9  2.054    2.074   -1.701     
Code
grouping_based_on_score = score_groups(as.data.frame(df_r3))
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = grouping_based_on_score) 
       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
Code
dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_r3, dif.sex)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.sex) 
       Item     Var   gamma     se pvalue padj.BH sig   lower   upper
1 MADRS.S_1 dif.sex  0.0377 0.0315 0.2313  0.6939     -0.0240  0.0995
2 MADRS.S_7 dif.sex  0.0597 0.0322 0.0637  0.1911     -0.0034  0.1229
3 MADRS.S_9 dif.sex -0.0840 0.0319 0.0084  0.0252   * -0.1465 -0.0215
Code
dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_r3, dif.age)

Item 3 4 5 Mean location StDev MaxDiff
MADRS.S_1 -0.058 -0.169 0.007 -0.073 0.089 0.176
MADRS.S_7 -0.243 -0.148 -0.304 -0.232 0.079 0.156
MADRS.S_9 0.301 0.317 0.297 0.305 0.011 0.020
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.age) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.1192 0.0181 0.0000  0.0000  *** -0.1546 -0.0838
2 MADRS.S_7 dif.age  0.0986 0.0184 0.0000  0.0000  ***  0.0625  0.1348
3 MADRS.S_9 dif.age  0.0290 0.0189 0.1259  0.3778      -0.0081  0.0661
Code
dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_r3, dif.edu)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 -0.086 -0.042 -0.064 0.031 0.044
MADRS.S_7 -0.254 -0.223 -0.239 0.022 0.031
MADRS.S_9 0.340 0.265 0.303 0.054 0.076
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.edu) 
       Item     Var   gamma     se pvalue padj.BH sig   lower  upper
1 MADRS.S_1 dif.edu -0.0443 0.0280 0.1137  0.3412     -0.0993 0.0106
2 MADRS.S_7 dif.edu  0.0240 0.0287 0.4028  1.0000     -0.0322 0.0803
3 MADRS.S_9 dif.edu  0.0120 0.0288 0.6775  1.0000     -0.0444 0.0683
Code
dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_r3, dif.yr)
[1] "No statistically significant DIF found."
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.yr) 
       Item    Var   gamma     se pvalue padj.BH sig   lower   upper
1 MADRS.S_1 dif.yr  0.0438 0.0192 0.0222  0.0667   .  0.0063  0.0814
2 MADRS.S_7 dif.yr  0.0233 0.0195 0.2328  0.6985     -0.0150  0.0615
3 MADRS.S_9 dif.yr -0.0674 0.0194 0.0005  0.0015  ** -0.1054 -0.0294

14.1 Analysis part 7 1 of 2

This was removing item 6

15 Rasch analysis 7 part 2 of 2

Removing item 7

Code
df_r3 <- df %>% filter(Time=='WEEK02') %>% select(!Time) %>% 
  select(!c(MADRS.S_4,MADRS.S_3,MADRS.S_2,MADRS.S_8,MADRS.S_5,MADRS.S_7))
Code
RIlistItemsMargin(df_r3, fontsize = 12)
itemnr item
MADRS.S_1 Sadness
MADRS.S_6 Lassitude
MADRS.S_9 Suicidal thoughts
Code
RIitemfit(df_r3, cutoff = "Smith98")
Item InfitMSQ Location 1 ± 2 / √n
MADRS.S_1 0.857 0.13 [0.973, 1.027]
MADRS.S_6 1.137 -0.54 [0.973, 1.027]
MADRS.S_9 1.031 0.50 [0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).
Code
simfit1 <- RIgetfit(df_r3, iterations = 300, cpu = n_cores)  

RIitemfit(df_r3, simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
MADRS.S_1 0.857 [0.966, 1.032] 0.869 [0.962, 1.033] 0.109 0.093 0.13
MADRS.S_6 1.137 [0.958, 1.034] 1.131 [0.961, 1.036] 0.103 0.095 -0.54
MADRS.S_9 1.031 [0.96, 1.038] 1.017 [0.945, 1.038] no misfit no misfit 0.50
Note:
MSQ values based on conditional calculations (n = 5465 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RIgetfitPlot(simfit1, df_r3)

Code
ICCplot(as.data.frame(df_r3), 
        itemnumber = 1:3, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[c(1,6,9)])

[1] "Please press Zoom on the Plots window to see the plot"
Code
### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)
Code
Item Observed value Model expected value Absolute difference Adjusted p-value (BH) Statistical significance level Location Relative location
MADRS.S_1 0.76 0.71 0.05 0.000 *** 0.13 2.10
MADRS.S_6 0.68 0.71 0.03 0.001 ** -0.54 1.42
MADRS.S_9 0.71 0.70 0.01 0.466 0.50 2.46
Code
RIbootRestscore(df_r3,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600
Item Item-restscore result % of iterations Conditional MSQ infit Relative average item location
MADRS.S_1 no misfit 32.0 0.86 2.10
MADRS.S_1 overfit 68.0 0.86 2.10
MADRS.S_6 no misfit 97.6 1.14 1.42
MADRS.S_9 no misfit 100.0 1.03 2.46
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.
Code
RIpcmPCA(df_r3)
PCA of Rasch model residuals
Eigenvalues Proportion of variance
1.66 57.6%
1.34 42.1%
0.01 0.2%
Code
simcor1 <- RIgetResidCor(df_r3, iterations = 500, cpu = n_cores)
RIresidcorr(df_r3, cutoff = simcor1$p99)
MADRS.S_1 MADRS.S_6 MADRS.S_9
MADRS.S_1
MADRS.S_6 -0.45
MADRS.S_9 -0.26 -0.51
Note:
Relative cut-off value is -0.348, which is 0.06 above the average correlation (-0.408).
Correlations above the cut-off are highlighted in red text.
Code
Item 1 Item 2 Partial gamma SE Lower CI Upper CI Adjusted p-value (BH)
MADRS.S_1 MADRS.S_9 0.275 0.026 0.224 0.326 0
MADRS.S_1 MADRS.S_6 0.16 0.028 0.106 0.214 0
MADRS.S_9 MADRS.S_1 0.116 0.027 0.064 0.169 0
Code
RIloadLoc(df_r3)

Code
mirt(df_r3, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code
# for fewer items or a more magnified figure, use:
#RIitemCats(df)
Code
# increase fig-height above as needed, if you have many items
RItargeting(df_r3)

Code

Code
iarm::score_groups(as.data.frame(df_r3)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)
  score_group    n
1           1 2816
2           2 2649
Code
dif_plots <- df_r3 %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code
clr_tests(df_r3, model = "PCM")

Conditional Likelihood Ratio Tests:
        clr df pvalue sig 
overall 243 16 0       ***
Code
Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.780    0.819   -3.489  -- 
MADRS.S_6  1.252    1.222    2.763  +  
MADRS.S_9  0.593    0.585    0.735     

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1  2.515    2.480    2.509  +  
MADRS.S_6  2.894    2.921   -1.975     
MADRS.S_9  1.992    2.000   -0.595     
Code
grouping_based_on_score = score_groups(as.data.frame(df_r3))
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = grouping_based_on_score) 
       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
Code
dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_r3, dif.sex)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.101 0.102 0.101 0.001 0.001
MADRS.S_6 -0.513 -0.701 -0.607 0.133 0.188
MADRS.S_9 0.412 0.599 0.506 0.133 0.188
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.sex) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.sex  0.1089 0.0301 0.0003  0.0009  ***  0.0499  0.1680
2 MADRS.S_6 dif.sex -0.0674 0.0286 0.0186  0.0557    . -0.1234 -0.0113
3 MADRS.S_9 dif.sex -0.0253 0.0307 0.4085  1.0000      -0.0854  0.0347
Code
dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_r3, dif.age)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.067 0.121 0.094 0.039 0.055
MADRS.S_6 -0.626 -0.554 -0.590 0.051 0.072
MADRS.S_9 0.559 0.433 0.496 0.090 0.127
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.age) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.0434 0.0177 0.0143  0.0430    * -0.0782 -0.0087
2 MADRS.S_6 dif.age -0.0271 0.0166 0.1021  0.3062      -0.0596  0.0054
3 MADRS.S_9 dif.age  0.0885 0.0179 0.0000  0.0000  ***  0.0534  0.1236
Code
dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_r3, dif.edu)

Item 2 3 Mean location StDev MaxDiff
MADRS.S_1 0.112 0.083 0.098 0.021 0.030
MADRS.S_6 -0.642 -0.454 -0.548 0.133 0.189
MADRS.S_9 0.530 0.371 0.450 0.112 0.159
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.edu) 
       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.edu  0.0468 0.0270 0.0828  0.2484      -0.0061  0.0998
2 MADRS.S_6 dif.edu -0.1020 0.0253 0.0001  0.0002  *** -0.1515 -0.0524
3 MADRS.S_9 dif.edu  0.0981 0.0272 0.0003  0.0009  ***  0.0448  0.1514
Code
dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_r3, dif.yr)

Item 3 4 5 Mean location StDev MaxDiff
MADRS.S_1 0.207 -0.387 0.043 -0.046 0.307 0.594
MADRS.S_6 -0.643 -0.068 -0.568 -0.426 0.312 0.574
MADRS.S_9 0.436 0.455 0.525 0.472 0.047 0.089
Code
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.yr) 
       Item    Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.yr  0.0309 0.0183 0.0918  0.2755      -0.0050  0.0669
2 MADRS.S_6 dif.yr  0.0341 0.0175 0.0517  0.1551      -0.0002  0.0685
3 MADRS.S_9 dif.yr -0.0688 0.0186 0.0002  0.0006  *** -0.1052 -0.0324

15.1 Analysis part 7 2 of 2

This was removing item 7

15.2 Analysis part 7 conclusion

Removing item 6 resulted in a better fitting scale.

16 Analysis conclusion

Resulting MADRS-S reworked scale thus has item 1,7 and 9.

16.1 Resulting items

Code
df_r3 <- df %>% filter(Time=='WEEK02') %>% select(!Time) %>% 
  select(!c(MADRS.S_4,MADRS.S_3,MADRS.S_2,MADRS.S_8,MADRS.S_5,MADRS.S_6))

what_scale <- gsub('\\.','',items_to_use)
Code
RIlistItemsMargin(df_r3, fontsize = 12)
itemnr item
MADRS.S_1 Sadness
MADRS.S_7 Inability to feel
MADRS.S_9 Suicidal thoughts

16.2 Reliability

Code
RItif(df_r3)

Code
RItif(df_r3,samplePSI=TRUE)

16.3 Person parameter

Code
RIpfit(df_r3)

16.4 Misfit check

Code
pfit_u3poly <- PerFit::U3poly(matrix = df_r3, 
                      Ncat = 7, # make sure to input number of response categories, not thresholds
                      IRT.PModel = "PCM")
misfits <- PerFit::flagged.resp(pfit_u3poly) %>% 
  pluck("Scores") %>% 
  as.data.frame() %>% 
  pull(FlaggedID)

misfits2 <- RIpfit(df_r3, output = "rowid")
Code
RIitemfit(df_r3, simcut = simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
MADRS.S_1 0.943 [0.966, 1.032] 0.956 [0.962, 1.033] 0.023 0.006 -0.05
MADRS.S_7 0.971 [0.958, 1.034] 0.967 [0.961, 1.036] no misfit no misfit -0.22
MADRS.S_9 1.096 [0.96, 1.038] 1.086 [0.945, 1.038] 0.058 0.048 0.33
Note:
MSQ values based on conditional calculations (n = 5465 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RIitemfit(df_r3[-misfits,], simcut = simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
MADRS.S_1 0.953 [0.966, 1.032] 0.904 [0.962, 1.033] 0.013 0.058 0.12
MADRS.S_7 0.984 [0.958, 1.034] 0.862 [0.961, 1.036] no misfit 0.099 -0.42
MADRS.S_9 1.085 [0.96, 1.038] 1.028 [0.945, 1.038] 0.047 no misfit 0.38
Note:
MSQ values based on conditional calculations (n = 5006 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RItif(df_r3[-misfits,])

Code
RIitemfit(df_r3[-misfits2,], simcut = simfit1)
Item InfitMSQ Infit thresholds OutfitMSQ Outfit thresholds Infit diff Outfit diff Location
MADRS.S_1 0.958 [0.966, 1.032] 0.96 [0.962, 1.033] 0.008 0.002 -0.09
MADRS.S_7 0.982 [0.958, 1.034] 0.973 [0.961, 1.036] no misfit no misfit -0.20
MADRS.S_9 1.064 [0.96, 1.038] 1.047 [0.945, 1.038] 0.026 0.009 0.35
Note:
MSQ values based on conditional calculations (n = 5229 complete cases).
Simulation based thresholds from 300 simulated datasets.
Code
RItif(df_r3[-misfits2,])

16.5 Item parameters

Code
Threshold 1 Threshold 2 Threshold 3 Threshold 4 Threshold 5 Threshold 6 Item location
MADRS.S_1 -4.46 -2.08 -0.01 0.38 2.35 3.52 -0.05
MADRS.S_7 -5.08 -2.60 -0.14 1.28 2.38 2.82 -0.22
MADRS.S_9 -3.54 -1.93 0.96 2.11 4.03 NA 0.33
Note:
Item location is the average of the thresholds for each item.
Code
item_param_matrix <- RIitemparams(df_r3,output='dataframe') %>% as.matrix(.)
write.csv(item_param_matrix,file=paste0("./results/item_params_",what_scale,'.csv'),row.names=TRUE)

16.6 Ordinal sum to interval score

Code
RIscoreSE(df_r3)
Ordinal sum score Logit score Logit std.error
0 -6.470 0.789
1 -4.981 1.002
2 -4.021 0.951
3 -3.219 0.893
4 -2.494 0.867
5 -1.767 0.846
6 -1.012 0.811
7 -0.358 0.771
8 0.174 0.735
9 0.659 0.703
10 1.131 0.680
11 1.586 0.669
12 2.016 0.672
13 2.437 0.691
14 2.880 0.735
15 3.401 0.805
16 4.112 0.869
17 5.444 0.719
Code
RIscoreSE(df_r3,output='figure')

Code
sum_to_latent <- RIscoreSE(df_r3,output = 'dataframe')
write.csv(sum_to_latent,file=paste0('./results/ordinal_sum_to_latent_',what_scale,'.csv'))

16.7 Thetas

The ordinal sum to interval score contains the information to transform into the below thetas, but we plot them here for convinience (based on 3 items)

Code
est_thetas2 <- RIestThetas(df_r3, method = "WLE")
hist(est_thetas2$WLE, 
     col = "#009ca6", 
     main = "Histogram of person locations (thetas)", 
     breaks = 17)

17 Software used

Code
pkgs <- grateful::cite_packages(cite.tidyverse = TRUE, 
                      output = "table",
                      bib.file = "grateful-refs.bib",
                      include.RStudio = FALSE,
                      out.dir = getwd())
formattable(pkgs, table.attr = 'class=\"table table-striped\" style="font-size: 15px; font-family: Lato; width: 80%"')
Package Version Citation
base 4.4.1 @base
car 3.1.3 @car
doParallel 1.0.17 @doParallel
easyRasch 0.3.3 @easyRasch
eRm 1.0.6 @eRm2007a; @eRm2007b; @eRm2009c; @eRm2013d; @eRm2015e; @eRm2019f
formattable 0.2.1 @formattable
furrr 0.3.1 @furrr
ggrepel 0.9.6 @ggrepel
glue 1.8.0 @glue
gridExtra 2.3 @gridExtra
iarm 0.4.3 @iarm
janitor 2.2.0 @janitor
kableExtra 1.4.0 @kableExtra
knitr 1.49 @knitr2014; @knitr2015; @knitr2024
matrixStats 1.4.1 @matrixStats
mirt 1.43 @mirt
mokken 3.1.2 @mokken2007; @mokken2012
patchwork 1.3.0 @patchwork
PerFit 1.4.6 @PerFit
psych 2.4.6.26 @psych
psychotools 0.7.4 @psychotools2021; @psychotools2022; @psychotools2024
psychotree 0.16.1 @psychotree2010e; @psychotree2011a; @psychotree2015b; @psychotree2018c; @psychotree2018d
reshape 0.8.9 @reshape
rmarkdown 2.29 @rmarkdown2018; @rmarkdown2020; @rmarkdown2024
tidyverse 2.0.0 @tidyverse

18 Session info

Code
R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Stockholm
tzcode source: system (glibc)

attached base packages:
 [1] parallel  grid      stats4    stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] gridExtra_2.3     knitr_1.49        car_3.1-3         carData_3.0-5    
 [5] grateful_0.2.10   easyRasch_0.3.3   doParallel_1.0.17 iterators_1.0.14 
 [9] furrr_0.3.1       future_1.34.0     foreach_1.5.2     ggdist_3.3.2     
[13] janitor_2.2.0     iarm_0.4.3        hexbin_1.28.5     catR_3.17        
[17] glue_1.8.0        ggrepel_0.9.6     patchwork_1.3.0   reshape_0.8.9    
[21] matrixStats_1.4.1 psychotree_0.16-1 psychotools_0.7-4 partykit_1.2-23  
[25] mvtnorm_1.3-2     libcoin_1.0-10    psych_2.4.6.26    mirt_1.43        
[29] lattice_0.22-5    eRm_1.0-6         lubridate_1.9.4   forcats_1.0.0    
[33] stringr_1.5.1     dplyr_1.1.4       purrr_1.0.2       readr_2.1.5      
[37] tidyr_1.3.1       tibble_3.2.1      ggplot2_3.5.1     tidyverse_2.0.0  
[41] kableExtra_1.4.0  formattable_0.2.1

loaded via a namespace (and not attached):
  [1] splines_4.4.1        bitops_1.0-6         R.oo_1.27.0         
  [4] cellranger_1.1.0     deSolve_1.40         rpart_4.1.23        
  [7] lifecycle_1.0.4      rprojroot_2.0.4      globals_0.16.3      
 [10] MASS_7.3-61          backports_1.5.0      magrittr_2.0.3      
 [13] vcd_1.4-13           Hmisc_5.2-1          rmarkdown_2.29      
 [16] yaml_2.3.10          sm_2.2-6.0           sessioninfo_1.2.2   
 [19] pbapply_1.7-2        RColorBrewer_1.1-3   abind_1.4-8         
 [22] audio_0.1-11         expm_1.0-0           quadprog_1.5-8      
 [25] R.utils_2.12.3       RCurl_1.98-1.16      pracma_2.4.4        
 [28] nnet_7.3-19          listenv_0.9.1        testthat_3.2.2      
 [31] RPushbullet_0.3.4    vegan_2.6-8          parallelly_1.40.1   
 [34] svglite_2.1.3        permute_0.9-7        codetools_0.2-19    
 [37] xml2_1.3.6           tidyselect_1.2.1     farver_2.1.2        
 [40] base64enc_0.1-3      jsonlite_1.8.9       polycor_0.8-1       
 [43] ks_1.14.3            progressr_0.15.1     Formula_1.2-5       
 [46] survival_3.7-0       systemfonts_1.1.0    tools_4.4.1         
 [49] gnm_1.1-5            ragg_1.3.3           snow_0.4-4          
 [52] Rcpp_1.0.13-1        mnormt_2.1.1         admisc_0.36         
 [55] xfun_0.49            here_1.0.1           mgcv_1.9-1          
 [58] msm_1.8.2            distributional_0.5.0 ca_0.71.1           
 [61] withr_3.0.2          beepr_2.0            fastmap_1.2.0       
 [64] fansi_1.0.6          digest_0.6.37        timechange_0.3.0    
 [67] R6_2.5.1             textshaping_0.4.1    colorspace_2.1-1    
 [70] R.methodsS3_1.8.2    inum_1.0-5           utf8_1.2.4          
 [73] generics_0.1.3       renv_1.0.10          data.table_1.16.4   
 [76] SimDesign_2.17.1     htmlwidgets_1.6.4    pkgconfig_2.0.3     
 [79] gtable_0.3.6         lmtest_0.9-40        pcaPP_2.0-5         
 [82] brio_1.1.5           htmltools_0.5.8.1    scales_1.3.0        
 [85] snakecase_0.11.1     irtoys_0.2.2         rstudioapi_0.17.1   
 [88] tzdb_0.4.0           checkmate_2.3.2      nlme_3.1-165        
 [91] curl_6.0.1           zoo_1.8-12           KernSmooth_2.23-24  
 [94] relimp_1.0-5         vcdExtra_0.8-5       fda_6.2.0           
 [97] foreign_0.8-86       pillar_1.9.0         vctrs_0.6.5         
[100] Deriv_4.1.6          cluster_2.1.6        dcurver_0.9.2       
[103] GPArotation_2024.3-1 htmlTable_2.4.3      evaluate_1.0.1      
[106] cli_3.6.3            compiler_4.4.1       rlang_1.1.4         
[109] future.apply_1.11.3  labeling_0.4.3       mclust_6.1.1        
[112] fds_1.8              plyr_1.8.9           rainbow_3.8         
[115] stringi_1.8.4        hdrcde_3.4           viridisLite_0.4.2   
[118] munsell_0.5.1        Matrix_1.6-5         qvcalc_1.0.3        
[121] hms_1.1.3            ltm_1.2-0            PerFit_1.4.6        
[124] readxl_1.4.3        

19 References